ilSIlET  I C 


f f£  8 


3 


AIR  FORCE  OFFICE  OF  SCIERIIFIC  KEoEAi\Oxi  (AiSC) 
NOTICE  OF  TRANSMITTAL  TO  DDC 

" .’a  Aochnical  report  has  been  revIo“;^d  nmi  ip 
o;  ■ roved  for  public  roloaso  1A>V  AIr»  ii)0—l«i  (i'b)« 
Mstrlbutlou  is  unlimited. 

A.  D.  BLOSE 

Technical  Information  Officer 


UNCLASSIFIED  C 

" — ■ ^ \ 

security  cl  ASSIPIC  op  ’’'his  page  ''When  Dstm  Fnref^.^1' 

1 REeasrpOCUMENTATION  PAGE 


iZ.  GOVT  ACCESSION  NO. 


^ TiVi  f rmnd  Subtitle)  ^ 

_KAGNETOSPHERIC  MAGNETIC  FIELD  MODELING^  / 


READ  TR VC  "IONS 

BEFORE  C(  -.IPLETING  FORM 

3.  RECIP(ENT'«  ,<»TAL0C  NUM8ER 


S.  TYPE  OF  REPORT  & PERIOD  COVERED 


Interim 


P.”/oison  ‘ 
K.  A./^fitzerj 


6.  PERFORMING  ORG.  REPORT  NUMBER 


s.  contract  OR  grant  NUMSERI-sJ 


F44620-75-C-OO33 


9.  performing  ORG-M2ATION  NAME  AND  ADDRESS 

McDonnell  Douglas  Astronautics  Company  - ^ 

5301  Bolsa  Avenue  / 

Huntington  Beach,  CA  92647 

II.  controlling  OFFICE  NAME  AND  ADDRESS 

AFOSR/NP 

Bolling  AFB,  Bldg  410  - 

Wash  DC  20332 

14.  monitoring  agency  name  a ADORESSCi/  ditietwn  from  ConttoUing  OUice) 


10.  PROGRAM  element,  project.  TASK 
AREA^  4'OWK  UNIT  KU^SERS 


/ 9751/05  / /'>  * 

61102F 

12..  BEOGRP-&ATE 


f3.  MuMBER  op  RACfiflT^  . — R— 

97^ 

IS.  security  CLASSTfo/  - 


UNCLASSIFIED 


1S«.  OECLASSIPICATION  DOWNGRADING 
SCHEDULE 


[16.  DISTRIBUTION  STATEMENT  (ol  thta  Report) 


Approved  for  Public  Release;  Distribution  Unlimited 


I 17.  DISTRIBUTION  STATEMENT  Co/  tho  mbetrmct  entered  In  Block  20,  ii  dlHerent  from  Report) 


[u.  supplementary  notes 


19.  KEY  WORDS  fContirtue  on  reeeree  aide  ti  neceeeery  end  identify  by  block  number) 

magnetosphere,  magnetic  field,  charged  particles,  coordinate  systems, 
electric  field 


20.  ABSTRACT  fContinue  on  reverae  aida  if  neceaeery  end  identtfy  by  blocit  number) 

*A  quantitative  model  of  the  magnetospheric  magnetic  field  and  associated  pro- 
cedures for  accurately  cataloging  charged  particle  data  out  to  and  beyond 
geosynchronous  orbit  is  developed.  The  magnetic  field  model  incorporates  all 
major  magnetospheric  current  systems  and  is  valid  for  all  tilt  angles;  i.e., 
angles  of  incidence  of  the  solar  wind  on  the  dipole  axis.  Tne  model  accurately 
represents  the  total  magnetospheric  magnetic  field  for  conditions  of  low 
magnetic  activity  and  to  a geocentric  distance  of  15  earth  radii  or  to  the  mag- 
netopause. A new  (B,  I)  coordinate  system  is  develooed  to  nore  accurately 


DD 


form 
JAM  73 


EDITION  OF  I NOV  65  IS  OBSOLETE 


UJCLASSIFIED  / 

, WCUHITV  CkASUPICATION  OF  THIS  PACK  flWiMi  Om«  BhiIWmO 


Annual  Scientific  Report  for  Contract  F44G20-75-C-0033 
Sponsored  by  the  Air  Force  Office  of  Scientific  Research 


PRINCIPAL  INVESTIGATOR:  W.  P.  Olson 
CO-INVESTIGATOR:  K.  A.  Pfitzer 


MCDOmniELL.  DOUGLAS  ASTffO^AUTtCS  COMPA/Vy-  bVEST 

5301  flofsd  Avenue,  Huntington  Beach,  CA  9264 7 


TABLE  OF  CONTENTS 


% 


Page 


1.  INTRODUCTION . . . 1 

2.  r-IAGNETIC  FIELD  MODEL 3 

2.1  Model  Assembly 4 

2.1.1  The  procedure 4 

2.1.2  Tests  in  Construction  8 

2.2  Model  Features  and  Uses 14 

2.2.1  Magnetic  Field  Representations  14 

2.2.2  The  Ordering  of  Charged  Particle  Data 24 

2.3  Associated  Induced  Electric  Field,  Tj  28 

3.  THE  B,  L COORDINATE  SYSTEMS  FOR  PARTICLE  DATA  ORGANIZATION  38 

3.1  Generalization  to  the  Distorted  Field  41 

3.2  Tilt  Effect  in  B,  L 45 

4.  SUMMARY  AND  CONCLUSIONS 47 

5.  APPENDIX  A - COMPUTER  CODES 50 

5.1  Subroutine  Description  . . 50 

5.2  Sample  Output  Description  53 

5.3  Program  Listing 55 

5.4  Sample  Program  Output  79 

6.  APPENDIX  B - COORDINATE  TRANSFORMATIONS  AND  RELATED  ANALYTIC 

DERIVATIONS 87 

6.1  Coordinate  Transformation  87 

6.2  Time  Derivation  of  the  Tilt  Angle 95 

7.  REFERENCES 97 


* Section  1 


INTRODL'CTION 


There  are  many  reasons  for  DoD  interest  in  the  development  of  accurate  quantitative 
representations  of  various  environmental  parameters.  This  is  because  the  near 
earth  space  environment  influences  many  military  hardv.'are  systems.  In  many  cases, 
these  influences  can  be  designed  around  but  in  others  they  cannot.  In  such  cases, 
system  performance  can  vary  v/ith  the  behavior  of  this  environment. 


In  this  report  v;e  review  work  perform.ed  at  McDonnell  Douglas  Astronautics  Company 
during  the  past  year.  The  purpose  of  this  v.ork  was  to  develop  a quantitative 
model  of  the  magnetospheric  magnetic  field,  and  associated  procedures  for 
accurately  cataloging  charged  particle  data  out  to  and  beyond  geosynchronous 
orbit.  The  main  reason  for  the  development  of  this  model  and  procedure  v/as  to 
provide  a capability  for  mapping  charged  particle  fluxes  produced  by  natural  and 
manmade  sources  of  ionization.  In  particular,  it  had  been  found  that  attempts 
to  accurately  m.ap  charged  particle  fluxes  from  points  above  the  earth's  surface 
to  geosynchronous  orbit  v/ere  not  very  accurate.  One  of  the  main  sources  of  error 
was  felt  to  be  inaccuracies  in  the  magnetic  field  model.  Thus,  one  of  the  tests 
of  the  modelling  exercises  completed  here  will  be  comparisons  made  by  DoD  between 
m.odel  calculations  and  observed  particle  fluxes. 


It  is  intended  that  this  model  and  procedure  will  be  used  also  for  other  purposes. 
The  predecessor  magnetic  field  model  (Olson  and  Pfitzer,  1974)  has  been  distributed 
to  several  dozen  groups  throughout  the  world  and  has  been  used  for  a variety  of 
reasons  that  range  from  very  practical  (the  location  of  the  foot  of  the  field  line 
to  various  synchronous  orbiting  satellites)  to  other  more  basic  studies  that 
hopefully  will  shed  light  on  the  basic  structure  of  the  magnetosphere  and  its 


physical  processes.  It  has  been  the  intent  of  the  MDAC  group  to  produce  a set  of 
quantitative  environmental  models  that  can  be  used  together  to  accurately  specify 
and/or  predict  several  environmental  parameters.  In  this  role  the  magnetic  field 
model  described  here  will  be  used  as  an  input  in  conjunction  with  other  quantitative 
models  of  neutral  upper  atmospheric  density  and  ionospheric  electron  density. 

This  interaction  between  models  will  improve  the  density  models  since  the  corpus- 
cular energy  sources  which  influence  both  ionospheric  and  upper  atmospheric 
parameters  at  high  latitudes  are  constrained  to  move  along  magnetic  field  lines. 

A magnetic  field  model  capable  of  describing  variations  in  charged  particle 
parameters  associated  with  "tilt  angle"  changes  and  minor  variations  in  magnetic 
activity  should  be  helpful  in  locating  such  gross  environmental  features  as  the 
poleward  edge  of  the  mid  latitude  trough. 


The  reader  of  this  annual  report  is  cautioned  that  the  results  presented  here  are 
not  to  be  considered  the  final  description  of  the  magnetic  field  model  and  the  new 
B,  L coordinate  system.  Rather,  it  is  intended  that  the  magnetic  field  model 
together  with  its  tests  and  suggested  uses  and  the  B,  L coordinate  system  will 
both  be  described  in  papers  submitted  to  the  Journal  of  Geophysical  Research  or 
a similar  journal.  Also,  documentation  on  the  computer  programs  generated  during 
this  contract  will  be  described  after  further  interaction  with  DoD  personnel  in 
a separate  manual. 


Section  2 

THE  MAGNETIC  FIELD  MODEL 

In  our  previous  magnetic  field  model  (Olson  and  Pfitzer,  1974)  all  major 
magnetospheric  current  systems  acting  as  sources  to  the  total  magnetospheric 
magnetic  field  were  considered  together  for  the  first  time.  These  current  systems 
include  the  magnetopause  current  system,  formed  by  the  interaction  of  the  solar 
wind  particles  with  the  total  magnetospheric  magnetic  field,  the  tail  current 
system  produced  by  particles  flowing  across  the  tail  of  the  earth's  magnetosphere 
and  returning  on  its  boundary,  and  the  quiet  time  ring  current  system  produced 
by  the  motion  of  charged  particles  moving  within  the  magnetosphere.  The  magneto- 
pause current  is  formed  by  particles  which  interact  with  the  earth's  magnetic 
field  and  then  return  immediately  to  the  magnetosheath.  These  particles  interact 
with  the  current  system  only  for  fractions  of  seconds  as  they  are  deflected  by  the 
magnetospheric  magnetic  field.  The  particles  flowing  across  the  tail  of  the 
magnetosphere  which  form  the  tail  currents  remain  part  of  that  current  system  for  ^ 
as  long  as  several  days  as  they  make  their  journey  across  the  tail.  The  quiet  time 
ring  current  system,  however,  is  formed  by  particles  trapped  along  magnetic  field 
lines  in  the  inner  magnetosphere.  These  particles  may  contribute  to  the  currents 
for  as  long  as  weeks  or  even  months.  Although  the  magnetopause  currents  flow  only 
on  a surface  (on  the  magnetopause  itself)  both  the  tail  current  system  and  the 
ring  current  particles  are  distributed  over  a large  region  of  space  in  the  magneto- 
sphere. This  had  previously  posed  a severe  limitation  on  magnetic  field  models  ‘ 

since  typically  scalar  potentials  were  used  to  represent  the  magnetic  field.  A 
scalar  potential  representation  cannot  be  used  in  the  region  of  the  source  currents 
for  the  field.  Thus,  in  the  Olson  and  Pfitzer  model  (1974)  the  equivalent  of  a , 

vector  potential  was  used  for  the  field  representation.  This  allowed  for  the 

' ^ 


first  time  the  accurate  representation  of  the  currents  flowing  in  a distributed 
fashion  throughout  much  of  the  magnetosphere.  That  model,  however,  is  limited 
to  the  case  of  perpendicular  incidence  of  the  solar  wind  upon  the  earth  dipole 
axis.  In  order  to  extend  the  time  of  validity  of  the  magnetic  field  model  it  is 
required  that  it  be  valid  for  all  angles  of  incidence  of  the  solar  wind  on  the 
dipole  axis  ("tilt  angles")  and  for  all  values  of  magnetic  activity.  In  the 
current  modelling  exercise  the  first  of  these  constraints  has  been  removed.  The 
model  described  below  was  developed  such  that  it  represents  all  major  magneto- 
spheric  current  systems  and  is  valid  for  all  tilt  angles.  It  is  believed  that 
this  model  is  capable  of  accurately  representing  the  total  magnetospheric  magnetic, 
field  about  30-40  percent  of  the  time  - v/hen  magnetic  activity  is  low. 

2.1  Model  Assembly 

2.1.1  The  Procedure 

As  with  the  1974  model,  a quantitative  representation  of  the  various  important 
magnetospheric  current  systems  was  first  developed.  The  ring  and  tail  currents 
are  initially  represented  in  the  form  of  wire  loops.  In  Figures  1 and  2 projections 
of  the  wire  loop  onto  the  xz  plane  (in  solar  m.agnetospheric  coordinates)  are  shown 
for  tilt  angles  of  0 and  35  degrees.  It  is  seen  that  in  the  inner  m.agnetosphere 
there  is  a nest  of  three  sets  of  wire  loops  used  to  represent  the  ring  current. 

The  tail  currents  flow  almost  linearly  across  the  upper  and  lower  boundaries 
of  the  plasmasheet  in  the  distant  tail  while  close  to  the  earth  they  curve  around 
such  that  they  are  in  close  proximity  to  the  ring  current  loops.  The  return 
path  of  the  tail  currents  was  constructed  to  flow  approximately  on  the  boundary  of 
the  magnetopause  with  the  shape  similar  to  that  used  in  the  deterr-.ination  of  the 
magnetopause  currents.  A three-dimensional  view  of  these  current  systems  is  shown 
in  Figure  3. 


Figure  3.  Three-dimensional  view  of  the  current  system 


The  currents  flov/ing  on  the  magnetopause  and  the  associated  magnetic  field  were 
determined  using  procedures  similar  to  those  developed  by  Olson  (1969).  This 
procedure  permits  the  currents  to  be  determined  from  a known  magnetopause  shape. 
Instead  of  using  a self-consistent  shape  determined  by  the  pressure  balance 
condition,  a more  realistic  shape  based  on  observational  data  was  used.  The 
magnetopause  was  therefore  represented  as  being  more  flared  in  the  dawn  dusk 
and  tail  regions  and  as  having  a flatter  subsolar  area.  The  magnetic  fields 
from  these  current  systems  were  then  obtained  using  the  Biot-Savart  Law.  Various 
primary  tests  were  used  once  the  magnetic  fields  were  obtained  to  compare  the 
results  with  observations  of  the  magnetospheric  magnetic  field.  Differences 
v/ere  determined  and  used  to  calculate  changes  in  the  locations  and  strengths  of 
the  various  magnetospheric  current  systems.  It  is  noted  that  the  initial 
estimations  of  the  placement  of  the  wire  loops  was  almost  sufficient  and  little 
adjustment  was  necessary. 

2.1,2  Tests  in  Construction 

The  purpose  of  this  model  v/as  to  provide  an  accurate  representation  of  the  total 
magnetospheric  magnetic  field  out  to  and  beyond  geosynchronous  orbit.  For  this 
reason  we  have  determined  the  magnetic  field  only  out  to  a geocentric  distance  of 
15  earth  radii.  The  model  should  not  be  used  beyond  that  region  as  it  v/as  not 
constructed  with  tests  valid  outside  of  that  region. 

One  of  the  most  important  tests  of  any  m.agnetospheric  magnetic  field  model  is 
the  com.parison  of  model  results  with  the  observed  magnetic  field  along  the  earth 
sunline.  In  Figure  4 the  quantity  AB  is  shov/n  along  the  earth  sunline.  AB  is 
the  difference  between  the  total  observed  field  and  the  geomagnetic  dipole  field. 
It  is  that  contribution  to  the  total  field  produced  by  the  magnetospheric 


currents.  It  is  seen  that  the  observed  field  is  quite  structured  in  the  inner 
magnetosphere.  Older  magnetospheric  models  containing  only  the  magnetopause  field 
are  seen  to  not  be  capable  of  describing  this  structure.  Ihdeed  it  is  the  field 
from  the  quiet  time  ring  current  that  produces  most  of  this  structure.  It  is 
seen  that  the  new  model  does  an  excellent  job  of  representing  the  observed  AB 
variations  along  the  earth-sun  line. 

Note  that  the  one  discrepancy  betv/een  these  observations  and  our  new  model  is  in 
the  tail  region  with  X^j,^  ^ Vd  earth  radii  (SM  subscripts  refer  to  solar  magnetic 
coordinates).  If  contours  of  constant  AB  are  plotted  in  the  noon  midnight 
peridian  plane  (the  solar  magnetic  or  solar  magnetospheric  xz  plane)  this  apparent 
discrepancy  can  be  accounted  for.  The  model  results  are  shown  for  a tilt  angle 
of  0“  in  Figure  5.  These  results  can  be  compared  with  the  observed  AB  contours 
reported  by  Sugiura  and  his  colleagues  (1971)  and  reproduced  in  Figure  6.  It  is 
noted  that  there  is  excellent  agreement  in  the  placement  of  the  AB  = 0 contours. 

In  order  to  explain  the  apparent  discrepancy  mentioned  in  Figure  4 it  must  be 
pointed  out  that  the  AB  = 0 contour  closes  in  the  tail  along  the  X axis  at 
approximately  -14  earth  radii.  In  Figure  7 the  AB  = C contour  for  the  tilt  angle 
of  35“  is  seen  to  cross  the  X axis  at  approximately  -9  earth  radii.  This  suggests 
that  in  data  sets  which  contain  data  averaged  over  all  tilt  angles  the  AB  = 0 
contour  should  lie  somewhere  near  X = - 10  earth  radii  as  is  seen  to  be  the  case 
in  the  observations  reported  by  Sugiura,  et  al.,  1971  shown  in  Figure  4. 

Attempting  to  fit  the  AB  contours  reported  by  Sugiura,  et  al.,  led  to  an 
abnormally  large  field  in  the  near  earth  portion  of  the  tail  in  the  Olson-Pfitzer 
1974  model.  This  problem  has  been  brought  to  light  in  the  present  analysis  and 
has  beer  corrected  in  the  present  model,  which  for  the  0“  tilt  case  and  all 


Figure 


others  better  represents  the  average  magnetospheric  field  configuration.  Evidence 
of  the  abnormally  large  field  in  the  model  was  pointed  out  by  R.  J.  Walker  in 
his  review  of  magnetospheric  field  models  (1976).  This  is  most  readily  observed 
when  the  magnitude  of  the  total  magnetic  field  is  plotted  in  the  equatorial 
plane.  Contours  of  constant  total  magnetic  fields  are  shown  in  Figure  8.  They 
are  seen  to  vary  smoothly  in  the  tail  region  along  the  earth-sun  line.  Also, 
there  is  no  buildup  in  field  intensity  near  the  subsolar  region.  These  tests 
suggest  that  the  0°  tilt  versi on  of  the  present  model  is  now  very  close  to  the 


2.2  Model  Features  and  Uses 

Having  passed  the  preliminary  tests  discussed  above,  the  model  v/as  used  in  two 
general  ways  to  further  test  it  and  provide  information  on  various  magnetospheric 
I magnetic  field  and  particle  population  features. 

2.2.1  Magnetic  Field  Representations 

j Field  lines  in  the  noon  midnight  meridian  plane  are  shov/n  in  2 degree  intervals 

i in  Figures  9 and  10  for  0°  and  -35“  tilt.  Since  some  of  the  field  lines  on  the 

i 

1 nose  of  the  magnetosphere  penetrate  beyond  the  magnetopause,  there  is  some 

;j  difficulty  in  assigning  a value  for  the  last  closed  field  line.  If  the 

'M 

, ■ definition  is  the  last  line  to  return  to  the  other  hemisphere  within  the  magneto-  I 

^.  ] pause,  then  the  last  closed  field  line  has  a latitude  value  of  about  78“  v;hile  for  I 

j ‘ I 

^ J j a tilt  angle  of  -35“  the  latitude  of  the  last  closed  field  line  is  lowered  to  I 

approximately  76*.  I 


Figure  8.  Contours  of  constant  IB]  in  the  equatorial  plane. 


The  extent  of  the  field  lines. into  the  tail  of  the  magnetosphere  is  also  more 
pronounced  than  in  all  previous  models.  This  is  most  clearly  illustrated  in 
Figure  11,  taken  from  Walker's  review  paper  (1976).  Two  curves  have  been  added 
to  his  figure;  that  of  the  new  model  and  the  equatorial  crossing  of  a pure 
dipole  field.  The  figure  gives  the  magnetic  latitude  as  a function  of  the 
equatorial  intercept  of  the  corresponding  magnetic  field  line.  Thus,  for  example, 
the  equatorial  crossing  of  the  field  line  from  the  magnetic  latitude  of  about  75“ 
at  midnight  is  approximately  15  for  a dipole  field  while  for  our  new  model  the 
field  line  originating  at  the  earth’s  surface  at  approximately  70°  magnetic  latitude 
crosses  at  15  earth  radii  in  the  tail.  This  feature,  the  extension  of  field 
lines  far  into  the  tail  even  during  quiet  magnetic  conditions,  is  well  documented 
observationally. 


Another  test  of  any  magnetospheric  magnetic  field  model  is  provided  by  several 
of  the  barium  cloud  releases  that  have  been  made  in  the  magnetosphere.  In 
Figure  12  one  of  those  releases  is  shov;n  and  the  magnetic  field  line  illuminated 
by  the  barium  is  compared  with  field  lines  calculated  from  several  models  with 
the  field  location  and  initial  direction  given  at  one  end  of  that  portion  of  the 
line  suggested  by  the  location  of  the  barium  cloud.  It  is  seen  that  the  calcu- 
lations using  our  new  model  perform  significantly  better  than  the  older  models  and 
along  the  length  of  the  optical  track  are  well  within  the  error  bars  for  the 
position  of  field  line. 


One  of  the  long  awaited  uses  of  accurate  magnetospheric  magnetic  field  models 
that  include  tilt  angle  effects  is  the  variation  in  conjugate  point  locations 
that  must  occur  even  during  quiet  magnetic  conditions.  If  the  location  of  one 
foot  of  the  magnetic  field  line  is  held  constant  over  a 24-hour  period  the 
question  of  the  location  of  the  opposite  end  of  that  field  line  arises.  In 


latitude  versus  the  equatorial  intercept  of  the  corresponding 
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calculated  field  lines  with  barium  cloud  observations 
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Figure  13  the  variation  in  tire  conjugate  point  location  of  a particular  field 
line  is  shown  for  winter  and  summer  solstices  and  for  spring  equinox.  It  is 
seen  that  at  summer  solstice  and  near  spring  equinox,  the  conjugate  point 
variation  in  the  point  is  quite  small.  Hov.'ever,  near  winter  solstice  the  variation 
becomes  significantly  larger.  This  is  because  near  winter  solstice  this  particular 
point  is  near  the  last  closed  field  line  location  and  considerable  daily  variations 
occur  as  is  expected.  Typically  it  appears  that  the  conjugate  point  variations 
for  field  lines  at  auroral  latitudes  should  be  on  the  order  of  several  tens  of 
kilometers  over  a daily  period..  Hov.ever,  sometimes,  for  field  lines  moving  near 
the  last  closed  field  line  location,  the  variations  can  be  significantly  larger. 
These  are  results  predicted  by  the  m.odel.  One  additional  check  of  this  and  other 
magnetospheric  magnetic  field  models  might  be  provided  by  the  search  for  large 
conjugate  variation  during  quiet  magnetic  conditions. 


Another  important  use  of  accurate  magnetospheric  magnetic  field  models  is  the 
calculation  of  the  intercept  at  the  earth's  surface  of  field  lines  extending 
through  geosynchronous  orbit.  Such  computations  are  necessary  in  order  to 
coordinate  particle  and  field  work  done  simultaneously  with  ground  based  and 
synchronous  satellite  platform  sensors.  In  Figure  14  the  location  of  the  foot 
of  the  field  line  through  synchronous  orbit  is  shown  for  a dipole  only  field  and 
for  three  different  tilt  angle  values  as  a function  of  local  tim,e.  (Although 
the  actual  field  line  geometry  does  not  vary  in  this  manner,  such  a representation 
sheds  light  on  the  magnitude  of  the  daily  variations  in  field  line  intercepts 
with  the  earth's  surface  possible  for  various  synchronous  satellite  locations.) 
First  it  is  noted  that  the  foot  of  the  field  line  can  be  as  much  as  2“  lov/er 
than  that  of  a dipole  only  field.  Early  models  that  contained  only  the  contribu- 
tion from  the  magnetopause  currents  placed  the  foot  of  the  magnetic  field  line 
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at  a latitude  higher  than  that  of  the  dipole.  The  field,  line  intercept  from  some 
actual  synchronous  orbit  locations  are  shown  in  Figure  15.  It  is  seen  that  the 
largest  variation  is  near  summer  solstice  at  both  geosynchronous  longitudes.  It 
is  also  noted  that  the  daily  variation  can  be  as  large  as  about  150  kilometers  in 
latitude  and  a similar  extent  in  longitude. 

One  final  test  of  this  tilted  magnetospheric  magnetic  field  is  to  compare  model 
calculations  of  the  daily  variations  in  the  components  of  the  magnetic  field  at 
the  location  of  the  geosynchronous  satellite  ATS-1  with  the  daily  magnetic 
variations  observed  at  the  location  of  the  satellite.  Such  a comparison  is  shov/n 
in  Figure  16  for  the  months  of  January  and  June  where  the  values  of  the  field  at 
each  hour  were  averaged  for  the  5 most  quiet  days  during  the  month.  It  is  seen 
that  the  variations  agree  quite  well  in  both  form  and  magnitude.  Previous  models 
do  not  come  close  to  repeating  the  pattern  of  the  observed  variations. 

2.2.2  The  Ordering  of  Charged  Particle  Data 

In  addition  to  using  the  model  directly  to  compute  various  magnetospheric  magnetic 
field  features  of  interest,  it  can  be  used  together  with  the  proper  formalisms  to 
study  and  organize  both  low  energy  and  energetic  charged  particle  data.  The  main 
purpose  in  the  development  of  the  magnetic  field  model  was  to  use  it  together 
with  adiabatic  invariant  theory  to  provide  a more  accurate  means  for  organizing 
low  energy  charged  particle  data.  This  has  resulted  in  the  development  of  a new 
"coordinate"  system,  the  B,  L coordinate  system,  it  is  described  in  detail  in 
Section  3.  This  model  and  the  previous  Olson-Pfitzer  (1974)  model  have  both  been 
used  very  successfully  to  predict  and  interpret  solar  cosmic  ray  data.  The  1974 
model  was  the  first  magnetospheric  magnetic  field  model  capable  of  predicting 
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magnetic  cutoff  latitudes  for-  the  penetration  of  solar  cosmic  rays  close  to  those 
observed.  Prior  to.  the  introduction  of  this  model  many  theories  had  been  proposed 
to  explain  the  discrepancy  betv/een  the  observations  and  older  model  predictions. 

In  the  present  model  some  of  the  remaining  discrepancies  have  been  resolved. 

One  of  the  major  successes  of  the  Olson-Pfi tzer  (1974)  zero  tilt  model  was  the 
lov;ering  of  the  cosmic  ray  cutoffs  by  over  4 degrees  which  brought  experimental 
observations  and  theory  into  close  agreement.  The  zero  tilt  model  brought  the 
noon  and  midnight  cosmic  ray  cutoffs  into  complete  agreement.  However,  a problem 
remained  at  dawn  and  dusk  v/here  the  cutoffs  perdicted  by  the  model  were  3°-4°  too 
high.  It  was  hypothesized  that  this  discrepancy  was  the  result  of  an  excessively 
large  field  near  the  dav/n  dusk  flanks  of  the  magnetosphere. 


The  present  model  brings  the  magnetic  field  values  near  the  magnetosphere  boundaries 
into  better  agreement  with  observations,  and  is  successful  in  lowering  these 
Cutoffs.  Table  1 summarizes  the  cutoff  calculations.  The  table  shows  that  the 
observed  and  so  calculated  cosmic  ray  cutoffs  are  now  in  substantial  agreement  at 
all  local  times.  Figure  17  shows  a sample  cosmic  ray  trajectory  for  a 5 MeV 
proton.  The  proton  enters  the  daylit  dav/n  side  of  the  magnetosphere  and  enters 
the  polar  cap  at  a latitude  of  68“.  A survey  of  cosmic  ray  access  using  this  n.odel 
has  shown  that  all  particles  v/hich  enter  the  polar  cap  within  3°-4°  of  the  cutoff 
enter  the  magnetosphere  through  the  daylit  dawn  side.  The  particles  which  enter 
the  center,  or  high  latitude  polar  cap,  either  penetrate  directly  through  the 
weak  field  region  of  the  dayside  cusps  or  arrive  via  the  lobes  of  the  tail  field. 

The  data  in  Table  1 indicate  a very  smiall  tilt  dependence  for  the  cutoff  latitudes. 

The  tilt  dependence  is  less  than  a degree  at  dav/n  and  dusk,  a degree  or  two  at 
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midnight  and  possibly  as  much  as  3®  at  noon.  These  predicted  effects  are  much 
smaller  than  the  observed  storm  time  variation  of  the  cutoffs.  Thus,  to  study 
the  tilt  dependence  in  the  experimental  data,  the  storm  time  dependence  must 
first  be  removed.  Masley  (1975)  has  gathered  a large  set  (over  500  observations) 
of  cutoff  values  using  data  from  the  OGO-6  satellite.  Here,  orthonormal  least 
squares  fitting  techniques  v/ere  used  to  fit  a parabolic  surface  in  tilt  and  Kp 
through  the  data  set.  The  function  which  gave  the  smallest  rms  errors  is 

CUTOFF  = a-j  + 82^  + + 

where  u is  the  tilt  angle 

Kp  is  the  magnetic  index  and  the 

a^.  are  the  least  squares  determined  coefficients 

Figures  18  and  19  show  the  comparison  between  theory  and  observation  at  Kp  = 1. 

The  small  dots  indicate  the  best  fit  line  v.'hen  K„  = 1,  the  large  dots  are  the 

p 

data  for  0 < Kp  < 2 and  the  open  circles  are  the  tilt  dependent  cutoff  calculations. 
There  is  good  agreement  between  theory  and  experiment.  The  scatter  of  points  for 
the  noon  data  set  is  much  larger  indicating  that  the  Kp  dependence  does  not 
adequately  remove  the  dynamic  variations.  Future  studies  must  use  other  parameters 
such  as  to  remove  this  dependence. 

2 . 3 The  Associated  Induced  Electric  Field,  Tj 

Any  time-varying  magnetic  field  has  associated  with  it  an  induced  electric  field 
as  given  by  Maxwell's  equations.  It  is  of  considerable  interest  to  determine 
whether  or  not  the  induced  electric  field  associated  with  the  time-varying  magnetic 
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Figure  18,  Tilt  dependence  of  the  cosmic  ray  cutoffs  for  the  midnight  meridian  in  the  northern  polar  cap. 


Figure  19.  Tilt  dependence  of  the  cosmic  ray  cutoffs  for  the  noon  meridian  in  the  northern  polar  cap 


field  model  described  above  is  of  importance  for  the  motions  of  low  energy  charged 
particles  in  the  magnetosphere.  Although  several  sources  of  electric  fields  in 
the  magnetosphere  have  been  proposed,  this  is  the  first  one  that  can  be  associated 
with  an  accurate  quantitative  model.  This  particular  induced  electric  field  is 
also  of  interest  because  its  periodicity  is  so  large.  Other  induced  electric 
fields  associated  with  magnetic  storms  and  magnetospheric  substorms  have  been 
discussed  but  their  time  constants  typically  range  from  several  minutes  to  a 
fraction  of  an  hour. 


It  is  therefore  with  great  interest  that  we  report  our  preliminary  findings  of 
I this  induced  electric  field.  Its  magnitude  and  the  strength  of  its  components 

I are  shown  over  a 24-hour  period  at  the  location  of  the  synchronous  satellite 

j ATS-1  at  winter  solstice  in  Figure  20.  The  magnitude  of  this  electric  field  is 

found  to  be  a large  fraction  of  a millivolt  per  meter.  Fields  of  this  size  are 

‘ ! 

regularly  assumed  to  be  of  importance  to  the  plasma  that  fills  the  m.agnetosphere. 
Contours  of  constant  induced  electric  field  in  the  equatorial  plane  are  shown  in 
Figure  21  at  winter  solstice  at  universal  time  = 0.  The  field  intensity  peaks 

i 

somev/here  in  the  mid  magnetosphere  and  then  falls  off  again  toward  the  m.agneto- 
..  pause. 

I 

. I 

^ The  procedure  for  determining  this  induced  electric  field  is  now  discussed. 

'j  The  induced  electric  field  is  most  readily  determined  from  the  magnetic  vector 

potential  where  IT  = D x and  = - |^.  Actually  Tj  is  determined  from  the 
; product  t",  = - (|^)  (1^)  • The  values  of  t were  determined  at  the  same  time 

i that  the  vector  magnetic  field  was  being  constructed.  Like  C,  is  given 

i f 

[ ! analytically  in  terms  of  the  tilt  angle,  y.  The  analytic  formula  for  the  time 

derivative  of  the  tilt  angle  is  given  in  Appendix  B. 
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INDUCED  ELECTRIC  FIELD  AT  ATS-1  AT  WINTER  SOLSTICE 


HOURS  (U.T.) 


CONTOURS  OF  CONSTANT  INDUCED  FIELD 
INTENSITY  AT  WINTER  SOLSTICE  (UT=0)  IN  THE  Z=0  PLANE 


Figure  21 


t 
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It  niust  be  pointed  out  that  these  results  are  preliminary  and  that  the  induced 
electric  field  Ej  is  by  itself  only  an  estimate  of  the  real  electric  field 
produced  in  the  magnetosphere  by  the  wobble  of  the  earth's  dipole  axis  with 
respect  to  the  solar  wind.  A secondary  electric  field  is  produced  by  charge 
separation  along  magnetic  field  lines  due  to  the  plasma  immersed  in  the  m,agnetic 
field.  This  plasma  tends  to  cancel  an  electric  field  along  the  direction  of  the 
magnetic  field  lines.  Thus,  the  approxim.ate  condition,  T • B"  = 0,  applies  where 
r is  total  magnetosphere  electric  field  produced  by  the  wobbling  dipole.  E"  is 
therefore  represented  as  having  two  sources  (1)  the  induced  field  which  is  given 
in  terms  of  the  magnetic  vector  potential,  and  a second  source  that  is  derivable 
from  a scalar  potential  and  is  produced  by  the  charge  separation  along  magnetic 


field  lines,  thus  t 


9A 

9t 


- V4>  where  is  the  electric  scalar  potential.  The 

9(J) 


f • ^ = 0 condition  then  implies  that  • B = V4>  • 


is  the  unit  vector  in  the  direction  of  B and  s is  an  elem.ent  of  length  along  f. 
Thus,  E||  is  known  everyv/here,  if  (J  is  known  over  some  boundary  it  can  then  be 
found  everywhere  in  the  magnetosphere.  A knowledge  of  (})  and  its  gradients 
together  with  the  values  of  the  induced  electric  field  suffice  to  actually 
determine  the  total  electric  field  in  the  m,agne  to  sphere. 


This  problem  v/as  a biproduct  of  the  present  study  and  not  investigated  in  detail. 
It  is  however  called  to  the  reader's  attention  in  this  report  and  it  is  felt  it 
is  an  imiportant  problem  basic  to  the  study  of  ragnetospheric  physics.  Its 
implications  may  be  of  considerable  importance.  First,  the  presence  of  such  a 
field  demands  that  magnetic  field  lines  can  no  longer  be  censid cred  as  equi- 
potential  surfaces.  Indeed,  this  work  has  shown  that  the  potential  along  the 
magnetic  field  lines  extend iji o_ to  geosynchronous  orbit  can  vary  by  as  much  as 
7 or  8 ki lovolts.  The  most  important  censequen ce_ o f such  a field,  hov/ever,  is 
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that  it  can  accelerate  charged- particles  in  the  n-.aonetosphore  locally.  This  is 
one  of  the  key  problems  in  magnetospheric  physics,  to  understand  how  plasmas  in 
the  magnetosphere  gain  and  lose  energy.  Such  an  electric  field  source  should  at 
least  accelerate  some  of  the  particles  some  of  the  time.  It  would  be  expected 
that  it  would  be  those  particles  that  drift  around  the  magnetosphere  with  a period 
close  to  1 day.  The  solar  wind  plasma  with  proton  energies  typically  on  the  order 
of  1 kilovolt  takes  approximately  1 day  to  rotate  around  the  magnetosphere  at  the 
location  of  geosynchronous  orbit.  The  consequences  of  this  induced  electric 
field  thus  appear  to  be  very  exciting. 


Section  3 


THE  B,L  COORDINATE  SYSTEMS  FOR  PARTICLE  DATA  ORGANIZATION 

In  order  to  properly  utilize  this  magnetic  field  model  to  organize  charged 
particle  data  it  is  necessary  that  the  behavior  of  the  particles  in  the  magnetic 
field  be  properly  understood.  Adiabatic  invariant  theory  has  been  used  together 
with  magnetic  field  models  to  develop  "coordinate  systems"  for  the  meaningful 
representation  of  the  particle  data.  Northrup  and  Teller  (1960)  discussed  the 
use  of  adiabatic  invariants  in  describing  the  motion  of  charged  particles  in  a 
magnetic  field.  These  invariants  can  be  used  to  describe  charged  particle 
behavior  when  the  change  in  the  magnetic  field  over  one  cyclotron  radius  is  very 
small,  a condition  that  is  satisfied  by  almost  all  charged  particles  trapped 
in  the  earth's  magnetic  field. 

The  first  adiabatic  invariant  states 


C = constant 


(3.1) 


i.e.,  the  component  of  energy  perpendicular  to  the  magnetic  field,  Wp,  divided 
by  the  total  magnetic  field,  B,  is  a constant.  If  a is  the  pitch  angle  of  a 
particle  (i.e.,  the  angle  that  the  velocity  vector  of  the  particle  makes  with 
the  magnetic  field  vector)  then  Wp  = W sin  a where  W is  the  total  energy  of 
the  particle.  Particles  trapped  in  the  earth's  field  bounce  back  and  forth 
along  field  lines  between  two  points  called  mirror  points.  The  mirror  points 
are  the  points  where  a * 90®  and  Wp  = W,  and  therefore  the  magnetic  field  at 
the  mirror  points,  is 


^mir 


C 


(3.2) 
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That  is  for  a particle  of  fixed  energy  the  field  at  the  mirror  points, 

is  a constant  and  may  be  considered  an  alternate  form  of  the  first  invariant. 

The  second  invariant  for  a particle,  the  integral  invariant,  is  defined  as 

' ' f <’-3) 

^r  (B  . ) 

' mir' 

where  ds  is  the  differential  path  length  along  the  line  of  force  connecting 
the  two  mirror  points  r and  r'  and  is  the  parallel  component  of  the  magnetic 
field.  As  adiabatic  particles  move  in  the  earth's  magnetic  field  they  bounce 
back  and  forth  between  their  mirror  points  and  slowly  drift  perpendicular  to 
the  magnetic  field  in  such  a way  that  B^.^  and  I remain  constant.  In  the. 
earth's  magnetic  field  the  locus  of  points  in  space  that  have  the  same 
and  I form  a ring  in  each  hemisphere  and  the  field  lines  which  connect  these 
rings  form  a surface  on  which  particles  mirroring  at  B^^.^  are  trapped.  This 
surface  is  referred  to  as  a drift  shell.  In  general  particles  that  mirror  at 
different  values  of  B along  a line  of  force  will  not  follow  the  same  drift 
shell  (shell  splitting).  Mcllwain  (1961),  however,  showed  that  v/hen  external 
magnetic  field  sources  are  neglected  (he  used  the  Jensen  and  Cain  512  term 
expansion  to  represent  the  earth's  main  field)  this  shell  splitting  effect  is 
quite  small.  Mcllwain  made  use  of  the  fact  that  all  particles  which  initially 
drift  through  a point  will  always  be  found  on  approximately  the  same  shell  and 
labeled  these  shells  with  a unique  number,  L,  where  L is  approximately  the 
radial  distance  in  earth  radii  to  the  shell  at  the  equator.  Therefore,  a 
coordinate  system  which  identifies  the  shell  (L)  on  which  a particle  moves  and 
identifies  where  on  the  shell  it  mirrors  (B^j,^^)  uniquely  identifies  the  particle. 


In  general  L,  the  shell  parameter,  has  the  definition  L = f(B,  I)  but  Mcllwain 
showed  that  for  a dipole  field 


where  M is  the  earth's  dipole  moment  I and 
representation  of  the  earth's  main  field, 
dipole  field  expression  at  numerous  points 
L in  the  real  field  is  then  obtained  using 
used  the  function 


(3.4) 

B are  calculated  using  the  full  series 
The  function  f is  evaluated  using  the 
in  space.  An  analytic  expression  for 
least  squares  techniques.  Mcllwain 


n=6 

I X" 

n=0  " 


(3.5) 


where  X = In 


and  a^  are  fitted  by  the  least  squares  procedures. 


The  current  version  of  Mcllwain's  (B,  L)  coordinate  system  considers  only  the 
internal  magnetic  field  sources  and  thus  is  valid  near  the  earth  (L  ? 3)  where 
the  external  field  sources  are  negligible.  At  larger  distances  (L  5 3)  one 
must  include  the  external  field  sources  in  order  to  adequately  describe  the 
magnetic  field.  By  including  the  external  field  sources  in  the  definition  of 
(B,  L)  the  validity  of  (B,  L)  can  readily  be  extended  to  about  L = 5.  For  L < 5 
the  asymmetries  of  the  field  are  small  and  the  approximation  that  particles  on 
a given  field  line  will  follow  the  same  drift  shell  is  valid.  Roederer  (1969) 
first  showed  that  at  large  geocentric  distances  (L  > 5,  where  the  field  dis- 
tortions become  large)  particles  of  different  pitch  angle  should  follow  different 
drift  shells  and  that  it  is  important  to  include  shell  splitting  effects  in 


order  to  understand  particle  behavior  in  this  region.  Pfitzer,  et  a1.,  (1969) 
using  data  from  the  geosynchronous  satellite  ATS-1  and  the  elliptic  orbit 
satellite  OGO-3,  first  -emonstrated  experimentally  the  importance  of  shell 
splitting.  It  was  shown  that  at  synchronous  orbit  particles  passing  through  a 
given  point  and  having  different  pitch  angles  will  follow  drift  shells  that  may 
separate  by  as  much  as  one  Rq.  (See  Figure  22) 

3.1  Generalization  to  the  Distorted  Field 

The  first  step  in  generalizing  the  (B,  L)  coordinate  system  is  to  modify  the 
magnetic  field  code  to  combine  an  accurate  expansion  for  the  internal  field 
sources  with  an  expansion  for  the  external  field  sources.  The  currently  avail- 
able (B,  L)  program  uses  the  1960  Jensen  and  Cain  expansion.  We  have  selected 
IGS  75  (Barraclough  et  al.  1975),  a more  up-to-date  and  accurate  version 
of  the  internal  field  code  and  have  combined  it  with  the  external  field  expansion 
discussed  in  the  Section  2.0.  The  internal  magnetic  field  code  is  in  spheri- 
cal geographic  coordinates  and  the  external  code  in  cartesian  solar  magnetic 
coordinates  and  therefore  several  coordinate  transformation  (see  Appendix  B) 
were  developed  in  order  to  combine  the  fields.  This  total  magnetic  field  is 
then  used  for  all  (B,  L)  calculations. 

Several  methods  of  extending  (B,  L)  to  the  distorted  field  case  are  possible. 

No  matter  how  distorted  the  field,  it  is  always  possible  to  map  directional 
particle  fluxes,  F(a),  by  either  using  the  adiabatic  invariants  I) 

directly  to  organize  the  data  or  by  defining  a new  L (L*)  which  is  a function 
of  the  pitch  angle  of  the  measured  particle  as  we^l  as  its  location.  This 
coordinate  system  would  most  accurately  organize  directional  particle  data. 
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FIGURE  22  -Equatorial  intersection  of  rlrift  shells  (using  the  Mead-Williams 
magnetic  field  model)  for  prticles  having  pitch  angles  of  a=65° 
and  a=90°  at  the  location  of  the  satellite  are  shown  when  the 
satellite  is  positioned  at  fc.6  Rr,  0°  latitude,  and  30°  local 
time;  R3  and  Dj  are  input  parameter  to  tne  Mead-Williams  model 
(from  Pfitzer,  et  al.,  1969). 


If,  however,  only  omnidirectional  data  with  unspecified  pitch  angle  distributions 
are  available,  very  little  can  be  done  to  improve  (B,  L)  beyond  the  inclusion 
of  the  external  field  sources  in  the  calculation  of  L.  However,  when  the  pitch 
angle  distribution  is  known,  then  (B,  L)  can  be  modified  to  represent  the  average 
shell  on  which  the  observed  particles  move.  In  this  study  we  have  used  this 
third  approach.  We  have  modified  L and  defined  an  T which  represents  the 
average  L value  when  the  particle  flux  is  known  to  be  isotropic. 

Mcllwain's  program  for  calculating  (B,  L)  evaluates  the  second  invariant  by 

integrating  from  the  location  of  the  satellite  (i.e.,  in  Equation  (3)  is 

set  to  the  value  of  B at  the  location  of  the  measurement)  to  the  conjugate 

point.  This  gives  a value  of  I which  is  identically  accurate  for  particles 

having  a pitch  angle  of  a = 90“  and  thus  represents  one  extreme  of  the  pitch 

angle  distribution.  All  other  particles  with  different  pitch  angles  will  drift 

on  shells  either  outside  or  inside  the  a = 90“  shell.  By  modifying  the  B^^^ 

used  to  determine  I to  represent  a particle  that  drifts  on  a shell  near  the 

center  of  the  drift  shell  distribution  we  can  define  a B„.„  and  then  determine 

mir 

the  average  second  invariant,  T,  for  the  isotropic  flux  case. 

A modified  (B,  L)  thus  replaces  the  (B,  L)  system  where  L can  be  obtained  from 
the  following  procedure: 

1)  Determine  A number  of  drift  shells  were  computed  for  isotropic 

particle  distribution  in  the  combined  external  plus  internal  field  and 
the  following  algorithm  was  empirically  developed  for 


®mir  ” ^tnir  *^mir  ^ 


where  is  the  minimum  B along  the  field  line. 


mir  mir 


B . = B , 

mir  mir 

where 


if  EL  < 3 


0.311653 

B_.^ 


B.  =C*B.  <3*B. 

mir  mir  — mir 

where 


if  EL  > 3 


C = 1.3  if  EL  > 5 
C = 1 + 0.15  (EL  - 3)  if  EL  < 5 
2)  Determine  T using  the  equation 


T . f ''7, . iui  y'"  ds 


fl  ! 


where  A and  A'  are  the  points  on  the  field  line  where  B(s)  = 

3)  Substitute  F . and  T into 
mir 


®mir* 


M * I In 


(3.8) 


where  M is  the  dipole  moment  and  a^  are  the  coefficients  of  Mcllwain's 
expansion. 


The  above  procedure  has  the  property  that  as  drift  shell  splitting  becomes  small 
r approaches  L.  When  drift  shell  splitting  becomes  important  L will  represent 
the  average  drift  shell  more  accurately  than  L.  The  procedure  is  only  slightly 
more  complicated  than  the  present  (B,  L)  procedure  and  the  only  difference  in 
computer  time  results  from  evaluating  I over  a slightly  longer  path  and  from 
having  a more  complicated  expansion  for  the  magnetic  field.  Since  U is  deter- 
mined for  an  isotropic  particle  distribution,  it  will  be  less  accurate  for 
other  pitch  angle  distributions.  Also,  since  T represents  an  average  drift 
shell,  exceedingly  strong  (e.g. , discontinuous)  radial  gradients  in  the  particle 
distribution  cannot  be  mapped  accurately  with  U.  For  the  highly  anistropic 
cases  and  for  strong  radial  gradients  only  directional  measurements  can  be 
mapped  and  T is  not  the  most  useful  approach,  instead  L*  which  is  a modification 
of  L that  takes  into  account  the  pitch  angle  of  the  measurement  must  be  used. 

3.2  Tilt  Effects  in  B,  U 

In  the  previous  paragraphs  we  have  discussed  our  extension  of  the  Mcllwain 
(B,  L)  coordinate  system  to  a highly  distorted  field.  Since  the  dipole  of  the 
earth  is  inclined  to  the  rotation  axis,  the  tilt  angle,  v,  changes  slowly. 

This  causes  the  magnetic  field  at  large  geocentric  distances  > to  vary 
pariodically  with  a frequency  that  is  slow  compared  to  the  time  it  takes  all 
but  the  lowest  energy  particles  to  drift  around  the  earth.  Therefore,  the 
particles  adiabatic  invariants  remain  constant  during  this  change,  and  at  the 
end  of  a day  when  a 24-hour  cycle  has  been  completed  the  particles  (if  no  other 
disturbances  have  occured)  will  be  in  their  initial  shell  positions.  Therefore 
even  though  the  magnetic  field  is  chancing  with  time,  I (for  a given  particle 
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changing  tilt. 


1 

i 

1 

j Since  the  magnetic  field  at  a given  location  varies  with  the  tilt,  u,  the 

1 

I value  of  U will  move  inward  and  outward  over  the  24-hour  period.  Thus,  L 

I correctly  labels  the  constant  drift  shells  on  which  the  particles  move  during 

the  adiabatic  field  changes. 


Appendix  A describes  the  computer  programs  developed  to  evaluate  L,  as  well  as 
all  of  the  programs  for  determining  the  magnetic  field  at  a given  point  in 
space.  The  user  is  given  the  option  of  calculating  L as  defined  initially  by 
Mcllwain,  calculating  L using  the  external  plus  internal  field  or  calucating  L 
using  the  external  plus  internal  field. 
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Section  4 

SUMMARY  AND  CONCLUSIONS 

A new  maqnetospheric  magnetic  field  model  has  been  developed  and  described. 

It  is  valid  throughout  the  inner  magnetosphere  and  we  believe  it  currently 
offers  the  best  description  of  the  total  magnetospheric  magnetic  field  during 
quiet  magnetic  conditions.  It  is  valid  for  all  angles  of  incidence  of  the  solar 
wind  on  the  earth's  magnetic  dipole  axis.  Its  region  of  validity  is  out  to 
a geocentric  distance  of  15  earth  radii  or  to  the  magnetopause,  whichever  is 
smaller.  The  model  is  semiempri ical . That  is,  it  is  based  on  both  the  several 
physical  principles  that  govern  our  current  understanding  of  the  magnetospheric 

magnetic  field,  and  on  several  data  sets  that  describe  quantitatively  the  struc-  , 

i 

ture  of  the  magnetospheric  magnetic  field.  The  quantitative  model  representation  * 

i 

was  done  by  fitting  both  the  spatial  and  the  tilt  angle  coordinates  simultaneously. 

The  resulting  four-dimensional  model  represents  on  average  the  magnetic  field 
values  input  to  the  fitting  routine  to  approximately  15  percent  of  their  values 
with  an  average  error  over  all  three  magnetic  field  components  of  0.27  gammas. 


The  model  has  been  tested  in  several  ways  and  found  to  offer  an  accurate 
description  of  the  field  in  the  inner  magnetosphere.  The  model  has  already 
been  us^d  to  predict  the  earth  intercept  of  field  lines  from  various  geosynchro- 
nous satellites.  This  information  is  of  vital  interest  to  satellite  experimental 
groups  wishing  to  coordinate  their  work  with  ground  based  and  balloon  experiments. 
The  primary  reason  for  developing  the  model  was  to  combine  it  with  adiabatic 
invariant  theory  to  offer  a more  accurate  means  for  organizing  charged  particle 
data  out  to  and  beyond  geosynchronous  orbit.  This  has  resulted  in  the  development 
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of  a new  "coordinate  system,"  the  (B,  L)  coordinates,  which  will  be  useful  for 
such  particle  mapping  exercises. 

Because  this  magnetic  field  model  is  in  some  sense  time-variant  (as  the  tilt 
angle  changes  the  magnetic  field  changes  through  the  day),  there  is  associated 
with  it  an  induced  electric  field.  Although  this  induced  electric  field  was 
anticipated  and  vector  potential  representation  of  the  magnetic  field  components 
were  obtained  in  order  to  easily  calculate  this  induced  electric  field,  it  was 
not  the  principle  purpose  of  the  contract  research  and  therefore  has  not  been 
studied  in  detail.  It  is  believed,  however,  that  this  electric  field  may  exert 
considerable  influence  on  the  magnetosoheric  plasma  and  its  study  may  shed 
light  on  many  important  problems  in  magnetospheric  physics. 

The  magnetic  field  model  has  also  been  used  together  with  Lorenz  force  calculations 
to  study  the  cutoff  in  magnetic  latitude  of  the  region  of  penetration  of  solar 
cosmic  ray  particles.  With  the  introduction  of  this  model  there  is  no  longer 
any  discrepancy  between  model  calculations  and  observations.  Prior  to  the  recent 
earlier  work  of  Olson  and  Pfitzer  (1974)  the  "cutoff  problem"  was  considered 
to  be  one  of  the  key  problems  in  cosmic  ray  physics  and  several  theories  were 
promulgated  to  explain  the  difference  in  the  observed  and  calculated  cutoff 
latitudes.  It  therefore  may  be  said  that  use  of  this  model  has  led  to  the 
solution  of  this  classical  cosmic  ray  physics  problem. 

Complete  listings  of  the  model  magnetic  field  and  the  new  (B,  U)  coordinate 
system  routines  are  given  in  Appendix  A.  It  is  again  called  to  the  reader's 
attention  that  this  document  is  meant  to  be  only  what  it  says,  an  annual 


report.  The  final  descriptions  of  the  maQ^ctic  field  model  and  its  testing  and 
uses  win  eventually  be  provided  by  journal  articles  which  will  be  submitted 
shortly.  Also,  the  listings  provided  in  the  c'pendices  should  be  considered 
provisional.  A final  authoritative  description  of  the  model  and  the  (B,  0 
coordinate  system  will  be  available  after  some  final  inputs  are  obtained  from 
Don  personnel. 
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Section  5 
APPENDIX  A - 
COMPUTER  CODES 

The  complete  computer  program  developed  under  this  contract  effort  consists  of 
a set  of  nine  separate  routines.  These  routines  may  be  combined  in  S''veral 
different  v/ays  depending  on  user  needs.  The  four  most  important  configurations 
are  shown  in  Table  A-1.  The  program  listing  which  is  included  in  this  appendix 
is  typical  for  configuration  1.  A test  deck  which  exercises  all  of  the  subrou- 
tines is  supplied  with  configuration  1.  This  test  program  can  be  used  to  verify 
correct  operation  of  the  code  on  the  user's  computer.  All  subroutines  have  been 
written  in  ANSI  standard  FORTRAN  IV. 


A definition  of  each  subroutine,  the  calling  sequences,  restriction-  and 
limitations  are  contained  in  the  COMMENT  cards  associated  with  each  routine. 

5 . 1 Subroutine  Description 
A brief  description  of  each  subroutine  follows: 

IHVTST  is  a test  program  designed  to  exercise  the  subroutines.  It  permits 
the  user  to  verify  the  correct  operation  of  the  program  on  his  computer. 

FLDINT  is  the  field  line  integration  routine  which  calculates  L or  T 
depending  on  the  option  selected.  It  requires  all  other  subroutines 
for  correct  operation. 

BMNEXT  determines  the  total  magnetic  field  at  a point  by  combining  an 

internal  magnetic  field  code  given  in  geographic  spherical  coordinates 
with  an  external  magnetic  field  code  given  in  cartesian  solar  magnetic 
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structure  for  determining  L 
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coordinates.  BHNEXT  requires  the  internal  field  routine  SPHRCM,  the 
external  field  routine  BXYZMU,  and  the  tilt  angle  routine  ANGLE. 

BMNEXT  is  a generalized  total  field  routine  and  can  provide  input  and 
output  in  either  spherical  or  cartesian  geographic  coordinates. 

BXYZMU  is  the  tilt  dependent  external  magnetic  field  subroutine.  It 
supplies  the  magnetic  field  contribution  from  the  external  current 
sources  (ring,  tail,  and  magnetopause)  as  a function  of  position  and 
the  tilt  of  the  dipole  axis.  Input  and  output  are  in  solar  magnetic 
coordinates.  This  routine  is  completely  self-contained  and  requires  no 
other  routines  for  its  operation.  Hov/ever,  subroutine  ANGLE  may  be  used 
to  determine  the  tilt  angle,  p,  as  a function  of  time. 

SPHRCM  is  a modification  of  J.  C.  Cain's  (USGS)  fast  SPHRC  routines.  They 
have  been  shortened  from  a fourteenth  order  expansion  to  a nineth 
order  expansion  and  combined  into  a single  subroutine.  A smooth  trun- 
cation scheme  has  been  added  such  that  for  large  distances  fewer 
terms  are  used.  The  terms  are  turned  off  gradually.  The  version  used 
for  this  study  contains  the  IGS  75  (Barraclough,  et  al.,  1975)  coefficients 
through  n = 9.  SPHRCM  may  be  used  independent  of  all  other  routines. 

ANGLE  is  a subroutine  which  calculates  the  position  of  the  sun  in  the 
celestial  sphere  and  determines  the  tilt  angle,  u.  as  well  as  the 
rotation  sines  and  cosines  for  transforming  from  geomagnetic  to  solar 
magnetic  coordinates.  ANGLE  is  self-contained  and  may  be  used  indepen- 
dent of  all  other  routines.  Derivations  of  these  angles  are  given  in 
Appendix  B. 

INTERP  is  a parabolic  interoolation  routine  which  fits  a parabola  through 

three  points  and  evaluates  the  function  at  a specified  location.  It  ‘ 

i 

is  an  independent  routine.  k 
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MAXFIX  determines  given  and  The  algorithm  for  is 

given  in  Section  3.2. 

CARMEL  is  a copy  of  C.  E.  Mcllwain's  (l.'niversity  of  California,  San  Diego) 

code  for  determining  L.  Given  B and  I,  CARMEL  calculates  L;  or 

mir 

given  F . and  T,  CARMEL  determines  U.  CARMEL  may  be  used  independent 
mi  r 

of  all  other  routines. 

5.2  Sample  Output  Description 

The  test  program  (INVTST)  produces  a sample  output  used  by  calling  the  field 
line  integration  routine  (FLDINT)  at  various  locations  and  times  (see  Section  5.4). 
L is  determined  for  three  separate  cases:  (1)  Internal  field  only  (L  INT), 

(2)  internal  plus  the  external  field  (L  INT+EXT),  and  (3)  L which  uses  internal 
plus  the  external  field  (L  AVE). 


* 

! 


The  definition  of  the  column  headings  for  the  sample  output  are  given  below: 

LAT  The  geographic  latitude  in  degrees. 

L0NG  The  geographic  east  longitude  in  degrees. 

R The  geocentric  radial  distances.  R = 1.0  is  6371.2  km. 

DAY0FYR  The  day  of  the  year  (1-366).  January  1 is  DAY0FYR  * 1. 

L The  value  of  L or  U. 

BMIN  The  minimum  value  of  the  magnetic  field  along  the  field  line. 

MAGL0NG  The  magnetic  longitude  in  degrees  measured  eastward  of  the  meridian 
through  the  dipole  and  the  geographic  pole  (-69®  W longitude). 
MAGLAT  The  magnetic  latitude  in  degrees.  +90®  magnetic  latitude  is  at 
11.7®  North  and  69®  West  geographic. 


CPTIME  Central  processor  execution  time  in  seconds  on  a CDC  6600  for  the 
determination  of  a single  value  of  L.  Because  of  the  uncertain- 
ties involved  in  measuring  such  short  times,  the  times  may  be  in 
error  by  as  mush  as  +0.002  second.  These  tines  are  useful  in 
informing  the  user  of  the  relative  computer  times  required  for  the 
three  L options. 


I 
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5.3  Program  Listing 


1 1 


12 

10 

1 


PROGRftM  INVTSTC  I NPUT , OUTPUT , T ftP E6  = OUTP UT  ) 

OinENSION  1AR(3) 

DATA  lAR/  lOH  L I NT  ,10H  L INT-t-EXT.lOH  L AVE  / 

LN=100 

DO  1 10=1,90,89 
DA=ID 

DO  I IU=l,19,6 

UT=IU-1 

LN=100 

DO  1 IL=1,31,30 

FLAT=IL-1 

DO  1 ILG=1, 181,90 

XL0NG=ILG-1 

DO  1 IR=2,8,3 

R=Ifl 

DO  1 1C=1,3 
CALL  SECOND(TA) 

CALL  FlDINT(FLAT,XL0NG,R,19  75.,DA,UT,  IC-2,  EL  , BH,  X MLON  G , X PIL  AT  ) 

CALL  SECOND(TB) 

TC=TB-TA 

LN=LN*1 

IF  (LN.LT.58 ) GO  TO  10 
LN=0 

PRINT  11 

FORHATC IHI, 15X,*LAT  LONG  R DAYOFVR  TI 

• PIE  L BniN  nAGLONG  (lAGLAT  CPTIdE*, 

*//  ) 

FORPIATt  AlO,  10F12 .5  ) 

PRINT  12,IAR(IC),FLAT,XL0NG,R,DA,UT,EL,Bri,  XMLONG  , XHL  AT , TC 

CONTINUE 

END 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 
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SUBROUTINE  ELDINT(  X L AT , X LONG,  R , VR  , 0 A Y , T I HE , J SUT  CH  , EL  , B Bn  1 N , X riLONG , 
*XHL AT  ) 

FLOINT  CALCULATES  THE  SECOND  ADIABATIC  INVARIANT  BY  INTEGRATING  ALONG 
HAGNETIC  LINES  OF  FORCE 

IT  CAN  CALCULATE  EITHER  THE  OLD  INTERAL  FIELD  ONLY  L,  L nOOIFIEO 
BY  THE  EXTERNAL  FIELD,  OR  AN  AVERAGE  L BASED  ON  ISOTROPIC  PARTICLES 
CALLING  SEQUENCE  - INPUT  PARAnETER 
XLAT  GEOCENTRIC  LATITUDE  INDEGREES 

XLONG  GEOCENTRIC  LONGITUDE  EAST  OF  GREENUHICH  IN  DEGREES 
R GEOCENTRIC  DISTANCE  FROn  THE  EARTHS  CENTER  IN  UNITS  OF 

EARTH  RADII,  HE.  RE=6371.2  Kn 

YR  THE  YEAR  - USED  BY  SOnE  nAIN  FIELD  ROUTINES  TO  SET  THE 

COEFFICIENTS  (E.G.  JULY  1 5 , 1 9 64= 1 9 6 M . 5 A ) 

NOTE***  YR  SHOULD  BE  CHANGED  ONLY  EVERY  FEU  DAYS  OR  nONTHS. 
NEU  FIELD  COEFFICIENTS  HUST  BE  COMPUTED  FOR  EVERY  CHANGE 
IN  YR,  THIS  COULD  CAUSE  A LARGE  INCREASE  IN  COMPUTER  TIME. 
THE  EARTHS  FIELD  CHANGES  ONLY  ABOUT  .001  GAUSS/VEAR  AT  THE 
EARTHS  SURFACE. 

DAY  THE  DAY  OF  YEAR  (1.-366.).  DAY  IS  USED  BY  THE  DIPOLE  TILT 
ROUTINE. 

DAY  MUST  BE  A UHOLE  NUMBER  AND  DAY  1 IS  JANUARY  1 
TIME  UNIVERSAL  TIME  IN  HOURS  (0.000-24.0000) 

JSUTCH  =-l  COMPUTE  L USING  INTERNAL  FIELD  ONLY 

= 0 COMPUTE  L USING  INTERNAL  + EXTERNAL  FIELD 

= ♦1  COMPUTE  AVERAGE  L USING  INTERNAL  EXTERNAL  FIELD 
OUTPUT  PARAMETERS 

EL  THE  L VALUE  DETERMINED  AS  REQUESTED  BY  JSUTCH 

BBMIN  THE  MINIMUM  VALUE  OF  B ALONG  THE  FIELD  LINE  IN  GAUSS 

XMLONG  THE  MAGNETIC  LONGITUDE  MEASURED  EAST  OF  THE  MERIDIAN 
CONNECTING  THE  TUO  POLES  (APP.EAST  OF  69U  GEOGRAPHIC) 

XMLAT  THE  MAGNETIC  LATITUDE 

****NOTE**** 

SINCE  THIS  ROUTINE  USES  AN  ACTUAL  MAGNETO SP HER  I C MAGNETIC  FIELD, 

THE  FIELD  LINES  ARE  NOT  ALL  CLOSED.  THUS  L IS  DEFINED  ONLY  IN  THE 
INNER  MAGNETOSPHERE.  AN  ATTEMPT  TO  CALCULATED  L OUTSIDE  OF  THIS 
REGION  WILL  set  L EQUAL  TO  100.  (EL=100.) 

BBMIN  WILL  BE  SET  TO  THE  LOCAL  VALUE  OF  8. 


COMMON/ BXYZ CM/ YEAR,  DAYYR,UT,XMLG,KODE,JSU,INITL,XMLT 
DIMENSION  X(3),XL(3),Q(3),QL(3),B(3),BL(3),SS(3),RR<3),BB(100,4), 
*S(100),SDS(100),B1NTL(3>,XX(3) 

DIMENSION  SX J(  100  ) 

DATA  P2 9, OPT/. 29289322, 1.70710678/ 

DATA  PICON/. 01745329252/ 

SET  UP  INITIAL  CONDITIONS 
YEARsYR 
UTsTlME 
DAYYRrOAY 
K00E=1 
JSU= JSUTCH 

COSINE=COS( XLAT*PICON  ) 

XX(  I )= R* C05I NE*COS( XLONG *P ICON) 

XX(2 )=R*COSINE*S!N( XLONG *P ICON) 

XX( 3 ) = R*S I N(  XL AT*P I CON  ) 

INITL=0. 


r 


nil. 


N=2 

CALL  BnNEXT( XX, B, BB( 2, 1 ) ) 

IF( XMLG.LT.O.  ) xnLG=XHLG*360. 

XnLAT=XMLT 

XnLONG=XnLG 

C EXIT  IF  OVER  THE  POLAR  CAP  OR  DISTANCE  IS  TOO  I ARGE  OR  FIELD  IS  TOO 
C UEAK. 

IF(  ABS(  XflLT  ).  GT.  75  . . OR  .R  . GT.  12.  .OR  . BB(  2,  1 >.  LT.0.00025  ) GO  TO  330 

IN1TL=1 

NBf1IN  = 2 

riAXFLG=0 

XJ  = 0 

ISU  = 0 

IPASS=0 

EL  = 0 

DO  10  1=1,3 
BINTLI I ) = B( I ) 

X(  I )=XX(  I ) 

10  0(1  1 = 0. 

S(  2 1 = 0. 

C SET  UP  ERROR  CONSTRAINTS 
ERRR= . 005 
ERR=ERRR/10 . 

SER  = SQRT((X(1  )**2*X(2)**2«-X(3)»»2)»ERR  ) 

FER=ERR**0 . 25 

DEL=-2.5*FER 

DS=SER 

BrnAX  = BB(  2,  1 ) 

BmN=B(1AX 

C STEP  ONCE  IN  THE  INCREASING  FIELD  STRENGTH  DIRECTION  AND  SET  STEP  TO 


c 

INTEGRATE  IN  DECREASING  FIELD 
IF( XMLT. GT. 0.  )G0  TO  30 

DIRECTION 

20 

DEL=-DEL 

DS=-DS 

30 

DO  35  1=1,3 

35 

XL(  I )=X( I )+DS*B(  I )/BB( 2, 1 ) 
CALL  BnNEXTt XL, BL, 8B( 1 , 1 ) ) 

IF( BB( 1 , 1 ). LT. BB( 2, 1 ) > GO 

TO  20 

C DETERHINE  CURVATURE  AND  OETERfllNE  STEP  SIZE 
S(  1 )=DS 
AO  CURV=0. 

CURVnN=ABS(  DEL  )/(1.5*(X(l  )**2*X(2)**2*X(3)**2)**(  ,75)*FER) 
DO  50  1=1,3 

5 0 CURV  = CURV  + ( (B(I)/BB(N, 1 )-BL(I)/BB(N-l,l  ))/DS>*»2 

CURV  = SQRT( CURV  ) 

CURV  = Af1AXl  ( CURV,  CURVMN  ) 

C OEL/CURV 

DS  = SIGN(AnlNI(ABS(DS),2.e),0S> 

C SAVE  LAST  STEP 
60  00  65  1=1,3 

Xl.(  I )=X(  I ) 

0L(  I ) = 0(  I ) 

65  BL(  I >=B(  I ) 

C BEGIN  GILLS  METHOD  OF  RUNGE  KUTTA  INTEGRATION 

66  P50S=.5*DS 
P290S=P29»DS 
0P70S=0P7»DS 


-Ti 


SOS(N)=DS 
DO  70  1=1,3 

SS(  r '=P5DS*B( t )/BB( N, t ) 

RH(  I ) = SS(  I )-Q(  I ) 

X(  I )=X(  I )*RR(  I ) 

70  0(  t »=Q(  I )*3 . *RR(  I )-SS< r ) 

CALL  BriNEXK  X,  B,  6B(  N,2  ) ) 

DO  71  1=1,3 

SS<  n = P2  9 0S*B(  I >/BB(W,2  ) 

RR(  I ) = SS(  I )-P29»0(  I ) 

X(  I )=X(  I )«RR(  I ) 

71  Q(  I ) = 0(  I )*3 . *RR<  I 1-SS<  I ) 

CALL  BriNEXTl  X,  B,  BB<  N,  3 > ) 

DO  72  1=1,3 

SS(  t ) = 0P7  0S*B( I )/B6(N,3  ) 

RR(  1 ) = SSf  I )-0P7*0( I ) 

X(  I )=X(  I )*RR(  I ) 

72  0( I )=Q(  I )*3 . »RR(  I )-SS(  I » 

CALL  BMNEX7(  X, B, BB( N,9  ) ) 

DO  73  1=1,3 

SS(  I ) = P5DS*B(  I )/BB(N,<*  ) 

RR(  I ) = ( SS(  r >-0(  I ) >/3. 

X(  I )=X(  I li'RRl  I ) 

73  0(  I ) = Q(  I )^3. *RH(  ! )-SS<  I ) 

N=N*1 

S( N ) = S( N-1  )*DS 

CALL  BHNEXK  X,  B,  BB<  N,  1 ) ) 

C IF  N IS  TOO  BIG  STOP  CALCULATION 
IF( N. GE. 100  ) GO  TO  300 

C IF  DISTANCE  TOO  LARGE  OR  FIELD  TOO  SHALL  EXIT 

IF( X( 1 >*»2*X( 2 »»*2*X<  3 >**2. GT. 199. . OR. BB( N, I ). LT. 0. 00025  )G0  TO  330 
IF(  ISU)  >200,  100,  150 
100  1F(  BB(  N,  1 ).  GE.  BMIN  )G0  '^0  110 

BnlN=BB(N,  1 ) 

NBMIN=N 

110  IF< BBC N, 1 ). LT. BHAX  ) GO  TO  90 
C IF  BHAX  IS  EXCEEDED  STOP  INTEGRATION 
IFI MAXFLG. ED . 0 ) GO  TO  190 
120  n=i 

125  CALL  INTERPI  BB(  N-2,  1 ),  5<  N-2  ),  BPIAX,SF,  1) 

N=N-n 

C CHECK  TO  SEE  THAT  LAST  STEP  IS  CLOSE  TO  BHAX 
IFI AB5( SF-S( N ) ). LT. SER  ) GO  TO  160 
0S= .97*1 SF-S<  N ) ) 

ISU=1 

130  IFin.EQ.O  ) GO  TO  60 

00  190  1=1,3 
XI  I ) = XL(  I ) 

01  I ) = QL(  1 ) 

190  B(  I ) = BL(  I ) 

GO  TO  66 

150  IFl  BBI  N,  1 ) . GT  . BflAX  ) GO  TO  120 
P1=0 

GO  TO  125 

160  IFI N-3  )1  00,  1 70, 200 
170  IFI IPA5S . NE. 0 ) GO  TO  200 
C TOO  FEU  POINTS  TO  FIND  L GET  ONE  nOHE 


ISui- 

DS  = 0. 5*( S( N-1  )-S( N ) ) 

GO  TO  130 

180  IF(  IPASS.EQ.O  ) GO  TO  185 
182  SX J( 3 ) = SX J< 2 ) 

SX J( 2 ) = SX J(  1 ) 

SX J( 1 ) = SX J<  H ) 

BB( 3, 1 )=BB( 2,1  ) 

BBC  2,  1 )=BB(  1 , 1 ) 

BBC  1 , 1 ) = bBC  1 ) 

N=3 

GO  TO  215 

C CALCULATE  L FROM  BfllN  IF  POINT  IS  NEAR  BniN  OR  IF  INTEGRATION  FAILS 
185  XJ=0. 

EL  = C0.311653/BBnlN)**Cl./3.  ) 

RETURN 

190  riAXFLG=l 

CALL  lNTERPCBBCNBmN-1,1  ),SCNBMIN-1  >,  DUnnV,  BBHl  N,  -1  ) 

CALL  riAXFI  XC  BnAX,  BBHIN  > 

GO  TO  no 

C 

C FIELO  LINE  INTEGRATION  COMPLETE  CALCULATE  INVARIANT 
200  XA=0 
QA  = 0 

SOSV  = SDS( 2 ) 

SXJC  2 ) = 0. 

NN=N-1 

DO  210  K=2,NN 
05  = S0SC  K ) 

P5DS=.5»0S 

P29DS=P29*DS 

0P7DS=0P7*DS 

FACT0R=0. 

IFCBBCK,!  l.LT.BMAX  )F ACTOR  = SQRTC  l.-B3CK,l  )/BMAX) 

SArP5DS*FACT0R 
R A=SA-Q A 
XA=XA*R A 
0A=QA+3 . »R A-SA 
FACT0R=0. 

IFCBBCK,2).LT.BMAX)FACT0R  = SQRTC1.-BBCK,2  )/BMAX  ) 

SA=P29DS*FACT0R 
R A=S A-P29»0 A 
X A=X A«R A 
OAzQA+S . *R A-SA 
FACT0R=0 . 

IFCBBCK,3).LT.B"lAX)FACT0R  = SQATCl.-BBCK,3)/BriAX) 

SA=0P70S»FACT0R 
RA=S A-0P7*QA 
XA=XA«R A 
QA=0 A«3 . *R A-S A 
FACT0R=0 . 

IFCBBCK,**  l.LT.BMAX  )FACTOB  = SQRTC  1 .-BBCK,M  )7BnAX) 

SA  = P5  OS*F ACTOR 
R Arc  s A-Q A )/3 . 

X ArX A*RA 
OArO A«3 . *R A-S A 
SX JC  K*1  )=XA 


i 


»i 


s! 


I ! 


1:' 


► 


210  CONTINUE 

C INTERPOLATE  TO  IdPROVE  INVARIANT  VALUE 

215  CALL  INTERPI BB< N-2, 1 >, SX J( N-2  ), BHAX, XXJ,  N) 

X J=X J + ABS<  XX J > 

C IF  2ND  PASS  OR  BHAX  WAS  NOT  CHANGED  STOP  INTEGRATION 
IF( IPASS.NE.O  ) GO  TO  230 

IF(  ABS(  BB(  2,  1 )-BPtflX  )/BriAX  . LT.  ERR  ) GO  TO  230 

IP ASS=1 

SX J( 1 ) = SXJ( 3 ) 

BB(  1 , 1 )=BB(  3,  I ) 

S(  1 )=S(  3 > 

CALL  1NTERP(BB(2,1  ),S(2),BnAX,SF,2> 

IF( ABSI SF  ). LT. SER  ) GO  TO  182 

DEL=-DEL 

DS=-SDSV 

DO  220  1=1,3 

X( 1 )=XX(  I ) 

Q(  I ) = 0 

220  B(  I )=B1NTL(  I ) 

N=2 
ISU  = 0 
GO  TO  60 

C CALL  nCILUJAINS  L EXPANSION 
230  CALL  CARHELI BHAX, X J, EL  ) 

RETURN 

300  WRITE! 6, 310  ) 

310  F0RnAT(2OH  FIELD  LINE  TOO  LONG) 

C IF  EL  CANNOT  BE  CALCULATED  SET  BBHIN  TO  LAST  VALID  VALUE  OF  B 
C AND  SET  EL  TO  100. 

33  0 BBnlN  = SQRT( BINTLI 1 )**2  + B I NTL ( 2 )**2+BINTL( 3 )**2  ) 

EL=100  . 

RETURN 

END 
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SUBROUTINE  BRNE X T(  X X , B , Bn AG  ) 

SUBROUTINE  FOR  HAGNETIC  FIELD  HAIN  PLUS  EXTERNAL 
VERSION  9/1/76 

DETERniNES  THE  VECTOR  HAGNETIC  FIELD  IN  GEOGRAPHIC  COORDINATES 
USING  A SPHERICAL  COORDINATE  EXPANSION  OF  THE  EARTHS  FIELO  AND  A 
CARTESIAN  COORDINATE  EXPANSION  OF  THE  BOUNDARY,  TAIL  AND  RING 
CURRENT  FIELDS  IN  SOLAR  HAGNETIC  COORDINATES 

INPUT  AND  OUTPUT  ARE  VIA  THE  ARGUNENT  LIST  AS  WELL  AS  THRU  A COnnON 
BLOCK  BXVZCn.  UT,  KOOE,  DAVYR,  XnLG,  XHLT,  JSU,  AND  INITL  ARE 
TRANSniTTED  THROUGH  COnnON. 

COnnON  BLOCK  ARGUnENTS 

KODE=l  SPECIFIES  THAT  INPUT  AND  OUTPUT  ARE  TO  BE  IN  CARTESIAN  COOROS 
KODE=a  SPECIFIES  THAT  INPUT  AND  OUTPUT  ARE  TO  BE  IN  SPHERICAL  COORDS 
UT  IS  THE  universal  TIHE  IN  HOURS  (0.-29.) 

DAYYR  IS  THE  DAY  OF  YEAR  IN  DAYS  --  IT  IS  ONLY  USED  IF  A TILTED 
HAGNETOSHER I C EXTERNAL  FIELD  IS  USED 

THIS  deck  IS  SUPPLIED  WITH  A TILT  DEPENDENT  VERSION  OF  THE 
EXTERNAL  FIELD. 

JSU  LESS  THAN  ZERO  — GET  INTERNAL  FIELD  ONLY 

IF  GREATER  THAN  ZERO  ''ET  INTERNAL  PLUS  EXTERNAL  nAGNETIC  FIELO 
INITL  IF  equal  to  ZERO  SET  UP  THE  VALUES  FOR  XnLT  AND  XnLG 

IF  NOT  EQUAL  TO  ZERO  DO  NOT  CHANGE  THE  VALUES  OF  XHLG  AND  XHLT 
XnLG  THE  nAGNETIC  LONGITUDE  IN  DEGREES. 

69  DEGREES  UEST  GEOGRAPHIC  IS  ZERO  nAGNETIC  LONGITUDE 
XnLT  THE  nAGNETIC  LATITUDE  IN  DEGREES. 

ARGUnENT  LIST 

XX  IS  A THREE  DinENSION  INPUT  ARRAY  CONTAINING  THE  POSITION  IN 
GEOGRAPHIC  COORDINATES 
IF  KODE  = 1 

XX(1)=X,  XX(2)=Y,  XX(3)=Z.  X,Y,Z  ARE  IN  EARTH  RADI!  Z IS  THE 
GEOGRAPHIC  N.  POLE  X IS  IN  THE  nERIDIAN  OF  GREENWICH  Y=X  CROSS  Z 
IF  KODE  = 2 

XX(1)  = R,  XX( 2 >=C0LATITUE  IN  DEGREES,  X X ( 3 ) = LONG  I T U DE  IN  DEGREES 
R IS  IN  EARTH  RADII,  LONGITUDE  IS  0-360  ♦ IS  EAST  OF  GREENUHICH 
B IS  A THREE  DinENSIONED  OUTPUT  ARRAY  CONTAINING  VECTOR  FIELD 
IF  KOOE  = 1 

B(1)=BX,  B(2)=BY,  B(3)=BZ  B IS  CARTISIAN  AND  IN  GAUSS 
IF  KODE  = 2 

B(  1 )=BR  (OUTWARD),  B(2)=BTHETA  (SOUTH),  B(3>=BPHI  (EAST)  IN  GAUSS 
BnAG  IS  THE  nAGNITIDUE  OF  THE  FIELD  IN  GAUSS 
FOR  FURTHER  INFORnATION  CALL  OR  URITE  K.A.  PFITZER  AT  nCDCNNELL 
DOUGLAS  IN  HUNTINGTON  BEACH  CALIFORNIA  (719)  896-3231 
DinENSION  X( 3 ), B( 3 ), XX( 3 ) 

COnnON/BX YZCn/ YEAR, DAVYR,UT,XnLG,KODE,JSU, I NITL, XHLT 
ConnON  /GCon/  st, Ct, sp, CP, aor, bt, bp, br, NnAx 

DATA  PICON /5 7. 29577951/, SIND, COSO/. 2027872959, .9792228106/, 

•S69, C69/. 9335809265, .35836799 95 /,UTLST,DAYLST/2*123956./ 

IF  EITHER  TinE  OF  DAY  OR  DAY  OF  YEAR  HAS  CHANGED  --  UPDATE  THE  TILT 
IF( UT. EQ . UTLST . AND. DAYYR . EQ . DAYLST  ) GO  TO  1 
UTLSTrUT 
DAYLST=DAYYR 

CALL  ANGLE  ( T I L T , SP S , CPS  ) 

IF<  KODE. GT.  1 ) GO  TO  3 

Cartesian  entry 


X(  1 ) = XX(  1 ) 

X( 2 )=XX( 2 ) 

X(  3 )=XX(  3 ) 

R2  = X(  1 )»*2  + X( 2 )**Z 
R = 5QHT( X( 3 >*»2  + R2  > 

R2  = 5QRT<  R2  ) 

CT=X( 3 )/H 
ST=R2/R 
CP=X(  1 )/R2 
SP  = X( 2 )/R2 
GO  TO  5 

C SPERTCAL  ENTRV 
3 R=XX(  1 ) 

THETAG=XX( 2 )/PI CON 
PHIG=XX( 3 )/PICON 
CT=COS( THETAG  ) 

ST  = S I N( THETAG  ) 

CP=COS<  PHIG  ) 

SP=SIN(PHIG) 

X(  1 ) = R*5T*CP 
X( 2 ) = R»ST*SP 
X( 3 ) = R*CT 

C ROTATE  THE  X-AXIS  INTO  THE  DIPOLE  LONG  -69 

5 BX=0. 

BVrO. 

BZ  = 0. 

IF( JSW . LT . 0 . AND.  I NITL . NE . 0 ) GO  TO  9 
7 XPP  = X(  1 )»C69-X( 2 )*S69 

VPP=  X(  1 )*S69+X( 2 )»C69 

C ROTATE  THE  Z-AXIS  INTO  THE  DIPOLE  AXIS  (11.8) 

ZP=XPP*SIND+X(3)»C0SD 
XP=XPP*COSD-X( 3 )*SIND 
YPrYPP 

IF(  I Nl TL . NE. 0 ) GO  TO  6 
C IF  INITL  IS  ZERO  SET  UP  XHLG  AND  XnLT 
XriLG  = 0 . 

IF(  YP  . NE.  0 . . OR  . XP  . NE.  0 . ) XflL  G = P I CON  • AT  AN2  ( YP  , XP  ) 

XnLTr 90 . - ACOSI ZP/R  >*P I CON  ^ 

I F(  JSU  . LT  . 0 ) GO  TO  9 fi 

C ROTATE  X-AXIS  BACK  TOWARD  THE  SUN  APPROXIMATELY  -(180+UT-69)  a 

6 X(  1 )=XP*CPS-VP*SPS  I 

X( 2 )=XP*SPS  + YP»CPS  B 

X(3)=ZP  ^ 

C GET  THE  EXTERNAL  CONTRIBUTION  !; 

CALL  BXYZnuC X, B, TILT  ) I 

C ROTATE  THE  X-AXIS  TO  THE  DIPOLE  LONG  APPROXIMATELY  (180+UT-69)  I 

BMX=B(  1 )*CPS+B( 2 )*SPS  ' 

BMY  = -B(  1 )»SPS+B( 2 )*CPS 

BMZ=B(3)  J 

C ROTATE  THE  Z-AXIS  BACK  TO  GEOGRAPHIC  THE  Z-AXIS  (-11.8) 

BGMXr3MX»C0SD*BMZ*SIND 

BGMYrBMY 

BGMZ  = -BMX*S!ND4-BMZ*C0S0 

C ROTATE  X-AXIS  back  TOWARD  THE  GREENWICH  MERIDIAN 
BX=BGMX»C69+BGMY»S69 
BY=-BGMX*S69*BGMY*C69 
BZ=BGMZ 


r 


i 


- f ' 


C GET  MAIN  FIELD 
9 CONTINUE 
AQRr 1 . /R 

CALL  SPHRCn 

C COMBINE  MAIN  FIELD  AND  EXTERNAL  FIELD 
I F( KOOE . GT  . 1 » GO  TO  10 

C CONVERT  COMPONENTS  TO  CARTESIAN  COORDS  (X  TOWARD  GREENWICH  MERIDIAN) 
B(  1 )=(  BXi-CP»(  ST»BR  + CT*BT  )-SP»BP  )*0 .00001 
B( 2 ) = ( BV*SP»( ST*BR  + CT»BT  )+CP»BP  )*0. 00001 
B(  3 ) = ( BZ  + CT*BR-ST »BT  )*0 . 00001 
GO  TO  20 

C CONVERT  TO  SPHERICAL  COORDINATES 
10  B(1  ) = (BR->-(BX*CP  + BY*SP  )*ST+6Z»CT)*0. 00001 

B(  2 )=(  BT*(  BX»CP  + BY*SP  )»CT-BZ*ST  )*0. 00001 
B(  3 )=(  BP  + BY»CP'BX»SP  )»0. 00001 
20  BMAG  = S0RT(8(1  )**2+B(2)»*2  + B(3)*»2> 

RETURN 

END 
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SUBROUTINE  BX  VZriU(  XK  , BE,  T I LT  ) 

VERSION  U/01/T6 

THIS  SUBROUTINE  COMPUTES  THE  MAGNETIC  FIELD  FROM  SOURCES  EXTERNAL 
TO  THE  IARTH.  no  internal  FIELD  IS  INCLUDED  IN  THIS  ROUTINE. 

IT  INCLUDES  THE  FIELD  CnNTR I BUT  I ONS  FROM  THE  MAGNETOPAUSE  CURRENTS, 
AND  Currents  distributed  throughout  the  magnetosphere  ithe  tail  and 
RING  CURRENTS).  IT  IS  VALID  FOR  ALL  TILTS  OF  THE  EARTHIS  DIPOLE 
AXIS  and  is  valid  curing  quiet  magnetic  CONDITIONS. 

A GENERALIZED  PRTHONORMAL  LEAST  SQUARES  PROGRAM  UAS  USED  TO  FIT  THE 
COEFFICIENTS  OF  A POWER  SERIES  (INCLUDING  EXPONENTIAL  TERMS)  THROUGH 
FOURTH  ORDER  IN  SPACE  ANO  SECOND  ORDER  IN  TILT. 

THIS  EXPANSION  HAS  BEEN  OPTIMIZED  FOR  THE  NEAR  EARTH  REGION  AND 
IS  valid  to  15  EARTH  RADII.  OUTSIDE  OF  THIS  REGION  THE  FIELD 
DIVE°GE5  RAPIDLY  AND  A TEMPLATE  SETS  THE  FIELD  TO  ZERO. 

IN  ORDER  TO  IMPROVE  COMPUTATIONAL  SPEED  THE  FIELD  IS  SET  TO  ZERO 
BELOW  2 earth  RADII.  (IN  THIS  REGION  THE  MAIN  FIELD  DOMINATES  AND 
THE  VARIATIONS  EXPRESSED  BY  THIS  EXPANSION  ARE  NOT  SUFICCIENTLY 
accurate  TO  PREDICT  VARIATION  ON  THE  EARTHIS  SURFACE). 

calling  sequence***** 

INPUT  parameters 

XX  IS  A 3 dimensioned  array  giving  THE  POSITION  OF  THE  POINT 

UIHERE  THE  FIELD  IS  TO  BE  CALCULATED.  XX(1),  XX(2),  AND  XX(3) 
ARE  respectively  X,  Y,  AND  Z EXPRESSED  IN 

EARTH  RADDII  IN  SOLAR  MAGNETIC  COORDINATES.  Z IS  ALONG  THE 
NORTH  DIPOLE  AXIS,  X IS  PERPENDICULAR  TO  Z ANO  IN  THE  PLANE 
CONTAINING  THE  Z AXIS  AND  THE  SUN-EARTH  LINE,  Y IS  PERP.  TO  X 
AND  Z FORMING  A RIGHT  HANDED  COORDINATE  SYSTEM.  X IS  POSITIVE 
IN  THE  SOLAR  DIRECTION. 

TILT  IS  THE  TILT  OF  THE  DIPOLE  AXIS  IN  DEGREES.  IT  IS  THE 

COMPLIMENT  OF  THE  ANGLE  BETWEEN  THE  NORTH  DIPOLE  AXIS  AND 
THE  SOLAR  DIRECTION  (PSI).  TILT=90-PSI. 

OUTPUT  PARAMETERS 

BF  IS  THE  3 DIMENSIONED  ARRAY  CONTAINING  THE  COMPONENTS  OF  THE 

FIELD  IN  SOLAR  MAGNETIC  COORDINATES.  BF(1),  BF(2),  AND  BF(3> 

ARE  THE  BX,BY,  AND  BZ  COMPONENTS  IN  GAMMA. 


DIMENSION  BF(3),XX(3),AA(6R  ),BB(6H  ), 
*A(32),B(32),C(22),D(22),E(32),F(32), 
DATA  ( IT  AC  I ),  1=  1 , 32  ) /2,  1,2,  1,2,2,  1, 
*2,  1 , 2,  1 , 2, 1 , 2, 2, 2,  1 / 

DATA  ( ITBC  I ), 1 = 1 , 22  ) /2,  1 , 2, 1 , 2, 2,  1 , 
DATA  { 1TC(  I ), 1 = 1 , 32  ) /I , 2, 1 , 2,  1 , 1 , 2, 
*1,2, 1,2, 1,2, 1,1, 1,2/ 

DATA  (AACI),I=I,6H)/-2.26836E-02,-1. 
•-3.12195E-0R,  9.50629E-03,-2.91512E- 
*-9. 26978 E- 05,  1.62929E-08,-l .27599E- 

* 8 . 96680E-09, -5 . 55850E-05,  1.37909E- 

* 9 . 52869E-07 , -7 . 1 7669E-09 , 9.98627E 

* 1.6066RE-05, -2. 29957 E- 01,-1  93777E- 

* 1.6O65eE-O3,-9.O1190E-O7,-3.15O69E- 

*-I .80676E-07,-! . 12022E-03,  5.985b8E- 
*-l  .98  73  7 E-06,  9 . 839  7xE-  1 0 , - 7 99379E 
*-1.3l968E-05,-1.29729E-09,  1.92930E- 

* 1 . 38  1 86E-05 , -2 . 8 I 599E-08 , 7.96386  E- 

*-1.81802E-07,-1.00278E-03,  1.98792F- 

*-2.96079E-06,-3.95831E-10,  1,02p90E- 


CCC99  ),DD(99  ),EE(69  ),FF(69  >, 
TT(9  ),!TA(32  ),ITBC22  ),ITCC32  > 
2, 1,2, 1,2, 1,2, 1,2, 2,  1,2, 2, 2,1, 

2,2,2,1,2,1,2,1,2,1,2,2,2,1,27 

1,2, 1,2, 1,2, 1,2, 1,1, 2, 1,1, 1,2, 


0 1 863E 
06,-1  . 
09  , 1 . 


9 . 
3. 
1 . 
2 . 


08 

05 
09 

06 

07, -5  , 

09  , 3 . 

08, -1  . 
06  , 2 . 
06,-1  , 
05,-1  . 


-09.  3 

5731 7E 
90732E 
9 1 8 15  E 
33662E 
09903E 
031 25E 
90009E 
829  72E 
91 769E 
69909E 
1 6925  E 
907  1 6E 


92986E 
03,  8. 
06,-1  . 

05  , 1 . 

10,-5  . 
03,-9, 
09  , 9 . 

06  , 5 . 

06  , 7 . 

09,-5. 
08,  2. 
05,  1 , 

08,-9. 


00, 

62856E-08, 
65983E-02, 
59296E-08, 
97620E-02, 
15606E-07, 
92887E-09, 
1 6509E-09  , 
9 1 7 3 7 E-09  , 
3037  1 E-OB , 
95099E-09  , 
1 7556E-08, 
00855E-05  , 


TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 

TOTAL 


• 1 . ’58 1 8E-07/ 

DATA  ( BB(  I ), I = I , 6H  )/  9.H7753E-02,  I . 9 5 9 8 1 E- 0 9 . - 1 . 8 2 9 3 3 E ♦ 0 0 , TOTAL 

» 5.59882E-09,  5 . 0 36 6 5 E- 0 3 , -2 . 0 7 6 9 8 E-06 , 1 . 1 095 9 E -0 1 , - 3 . 9 5 8 3 T E-05 , TOTAL 

»-9  . HOlTBE-OS  , 5 . 06969E-07  , -I  . 20  11  2E-03,  3.699HE-06,  1.99899E-01,  TOTAL 
*-7 . 99929E-05 , 2.96382E-09,  9 . 6 5 8 7 0 E-O 7 , - 9 . 5 9 8 8 1 E- 09 , 2.9369  7E-07,  TOTAL 

• 3.06520E-09  , 3.07836E-07  , 6.9830  1 E-03,  1 . 2 6 2 5 1 E-0 6 , - 7 . 0 9 5 9 8 E-0 3 , TOTAL 

•-1  .55596E-05  , 3 . 0 6 9 6 5 E + 0 0 , - 7 . 8 9 8 9 3 E- 0 5 , 9.95  195E-03,  3.7  1 92  1 E-06,  TOTAL 

•-1.52002E-01,  6.81988E-06,-8.55686E-05,-9.01230E-08,'3.71958E-09,  TOTAL 

• 1 . 309  76E-07  , -1  . 8297  1 E-0  1 , 1 . 5 1 3 9 0 E-0 5 , - 1 . 9 5 9 1 2 E-09 , - 2 , 2 2 7 7 8 E -0 7 , TOTAL 

• 5 . 992 78E-05 , -3 . 72758E-06 , -1 . 5 9932 E-03,  8.09921E-06,  5.38012E-01,  10TAL 

*-l  .931  82E-09,  1.50000E-09  , 5 . 8 8 0 2 0 E-0 7 , - 1 . 5 9 0 00 E- 02 , 1.60799E-06,  TOTAL 

• 3.17837E-09,  1 . 78959E-07, -8 . 93799E-03,  6.37599E-06,  1.27887E-03,  TOTAL 

•-2 , 95  878E-07,-l  .932  1 0E-01  , 6 . 9 1 2 3 3 E-0 6 , - 2 . 8 0 6 3 7E-0 9 , -2 . 5 7 0 7 3 E -0 7 , TOTAL 

• 5.7839  3E-05  , 9.52128E-1  0,  1 . 8 9 6 2 1 E -09 , - 9 . 8 9 9 1 1 E-C 8 , - 1 . 5 005 8 E-0 2 , TOTAL 

• 6.21 772E-06 / 

DATA  (CC(I  ),r=l,99)/-1.88177E-02,-1.92993E-06,-2.89069E-01,  TOTAL 

*-8 . 99939E-05 , -9 . 76380E-09 , -9 . 5 2998E-08 , 1.61086E-03,  3.18728E-07,  TOTAL 

• 1.29159E-06,  5.52259E-10,  3.95593E-05,  5.61209E-08,  1.38267E-03,  TOTAL 

• 5.79237E-07,  1.86989E-06  , 7.1  01  75  E-1  0,  1 . 9 5 29 3 E-0 7 , - 2 . 9 7 5 9 1 E- 1 0 , TOTAL 

*-2.93029E-03,-6.70000E-07,-2.30629E-02,-6.22193E-06,-2.90815E-05,  TOTAL 

• 2.01689E-08,  1.76721E-09,  3.78689E-08,  9.88996E-06,  7.33820E-09,  TOTAL 

• 7.321  26E-05  , 8.93986E-08,  8 , 3 29 9 9 E- 06  , - 6 . 1 1 7 08 E-08  , 1.7888  1 E-09,  TOTAL 

» 8.62589E-07,  3.93729E-06,  2 . 5 3 7 8 3 E-0 9 2 . 09 2 39 E-0 7 , 8.16691E-10,  TOTAL 

• 1.68075E-05,  7.62815E-09,  2.26026E-09,  3.66391E-08,  3.99637E-07,  TOTAL 

» 2 . 2553 1 E-1 0/ 

DATA  (00(11,1=1,99)/  2.50193E-03,  1.01200E-06,  3.23821E+00,  TOTAL 

• 1.08589E-05,-3.39I99E-05,-5.27052E-07,-9.96161E-C2,-1.95913E-09,  TOTAL 

*-9  . 236  1 9E-06 , 1 . 9315  3E-08  , -2 . 6299f E-09  , 1 . 05  I 38 E-0 7 , -2 . I 5 7 8 9 E-0 1 , TOTAL 

*-2.2071  7E-0  7,-2.65687E-05  , 1.26370E-08  , 5 . 8 8 9 1 7 E- 0 7 , - 1 . 1 3 6 5 8 E-0 8 , TOTAL 

• 1.69385E-03,  1 . 99263E-06 , -1  . 66095E-0 1 , -1  . 96096E-05  , 1.228  1 1 E-09,  TOTAL 

• 3.93922E-08,  9 . 6 6 76 0 E-0 5 , - 6 . 3 2 1 5 0 E-0 7 , -9 . 9 7 9 00 E-05 , -2 , 7 8 5 7 8 E-0 8 , TOTAL 

» 1.77366E-02,  2.05901E-07,-1.91756E-03,-9.99392E-07,-1.99988E-01,  TOTAL 
*-2 . 071 70E-06 , -5 . 90993E-05 , 1.59289E-08,  7.30919E-05,  3.38786E-08,  TOTAL 

*-l ,59537E-09,-1 .65509E-07,  1.90990E-02,  2.03238E-06,  1.01198E-09,  TOTAL 

• 5.20815E-08/ 

DATA  ( EE(  1 ),  1 = 1 , 69  1/-2 . 77929E  + 0 1 , -1  . 0 1 957E-03  , 9.2  1 936E-02,  TOTAL 

*-8 . 521  77E-06,  5.1  91  06E-0  1 , 8 . 2 8 8 8 1 E-0 5 , -5  . 5 96 5 1 E-09  , 1.1  6736E-07,  TOTAL 

•-2. 11206 E-0 3, -5. 35969 E- 07,  9. 91990 E-01,-1.33679E-05, -7. 18692 E-09,  TOTAL 

» 6.  1 7358E-08,-3.5  1990E-03,-5  .29070E-07,  1 . 8 89 9 3 E-0 6 , -6 . 6 06 9 6 E- 1 0 , TOTAL 

*-l . 39 708E-03,  1.02160E-07,  1.58219E-06,  2.05090E-10,  1.18039E+00,  TOTAL 

• 1.58903E-09,  1 . 86999E-02,-9 .96977E-06,  5.99869E-02,  9.99690E-06,  TOTAL 

•-1 . 18335E-09,  6 . 95689E-09, -2 . 73839E-09 , -9 . 1 7883E-08 , 2.79126E-02,  TOTAL 

J *- 1 . 02567E-05 , -1  . 25927E-09  , 3 . 0 7 1 9 3E- 0 8 , -5 . 3 1 8 2 6 E- 09  , - 2 . 9 8 9 7 6 E-0 8 , TOTAL 

, »-9 . 89899E-05 , 9.91980E-08,  3.85563E-01,  9.16966E-05,  6.79799E-09,  TOTAL 

•-2.08736  E-07,-3.92659E-03,-3.13957E-06,-6.31361E-06,-2.  .92981E-09,  TOTAL 
•-2 . 63883E-03, - 1 . 32235  E-07, -6 . 1 9906E-06 , 3.59339  E-09,  6.65986E-03,  TOTAL 

*-5. 81999 E- 06,-1  .88809 E-09,  3.62055E-08,-9,69380E-09,-2.21159E-07,  TOTAL 

' *-l  . 77996E-09,  9 . 955  60E-08  , -3 . 1 886  7E-09 , -3 . 1 7697E-0  7, - 1 . 058 1 5 E-05 , TOTAL 

• 2.22220E-09/ 

DATA  ( FF(  I >,  1 = 1 , 69  1/-5  .07092E*00  , 9 . 7 1 9 6 0 E-0 3 , - 3 . 79 85  1 E- 0 3 , TOTAL 

. *-3 . 67309E-06 , -6 . 02939E-0 1 , 1.08990E-09,  5.09287E-09,  5.62210E-07,  TOTAL 

• * 7.05  7 1 8E-02,  5 . 1 3 1 6 0 E-0 6 , - 2 . 8 5 5 7 1 E + 0 0 , -9 . 3 1 7 2 8 E-05  , 1.03I85E-03,  TOTAL 

i • 1.05332E-07,  1.09106E-02,  1.60799E-05,  9.18031E-05,  3.32759E-08,  TOTAL 

? * I.20113E-01,  I . 90986E-05 , -3 . 37993E-05 , 5.98390E-09,  9.10815E-02,  TOTAL 

j *-9 . 00608E-09,  3 . T5393E-03, -9 . 69939E-07 , -2 . 98561 E-02 , 1.31836E-09,  TOTAL 

: , »-2 . 6 7755E-09  , -7 . 60285 E-08 , 3 . 0 9 9 9 3 E- 0 3 , - 3 2 8 9 5 6 E -0 6 , 5.8  2 36  7E-01,  TOTAL 

I • 5 . 39996E-06 , -6 . 1 526 1 E-09 , 9.05316E-08,  1 . 1 35 9 6 E- 0 2 , -9 . 2 6 9 9 3 E- 0 6 , TOTAL 


^ r 


•-2  . T?00  re-Oa  , 5 . 2i;5  23E'08  , '2 .985  76E  + 00,  3.0/325E-05,  1.516458-03,  TOTAI. 

• 1.250988-06,  9.072138-02,  1.059698-05,  1.092328-09,  1.773818-08,  TOTAL 

• 1.927818-01,  2.157398-05,-1.657418-05,-1.886838-09,  2.448038-01,  TOTAL 

• 1.513168-05,-3.011578-09,  8.970068-08,  1.869718-02,-6.990798-06,  TOTAL 

• 9.131988-03,-2.380528-07,  1.285528-01,  6.925958-06,-8.367928-05,  TOTAL 

•-6 . 100218-08/ 

DATA  TlLTL/99./ 

C 

X = XX(  1 1 
V=XX{ 2 ) 

2=  XX( 3 ) 

V2= 

Z2=2»»2 
R2=X»*2+V2+22 
BX  = 0. 

BV  = 0 . 

BZ  = 0. 

C TEST  FOR  LOCATION  OF  INPUT  POSITION 

C IF  TH8  location  is  0UTSID8  THE  REGION  OF  SERIES  CONVERGENCE  SET 
C THE  HAGNErrC  FIELD  TO  ZERO  AND  PRINT  AN  ERROR  HESSAGE 
C0N=1 . 

I F( R2 . GT . 225 . > GO  TO  50 
IF( R2 . LT. 9 . )G0  TO  90 

C TAPER  THE  FIELD  LINEARLY  TO  ZERO  IF  INSIDE  R = 2.5. 

IF(R2.LT.6.25  1 CO N=C0N»<R2-9. 01/2.25 
I F( T 1 LTL . EQ . T I LT  IGO  TO  6 
TILTL=TILT 

C IF  THE  TILT  HAS  CHANGED  UPDATE  THE  CO EFF I C I ENTS . 

TT(  I 1=1 
TT( 2 ) = TILTL 
TT( 3 )=TILTL»*2 
TT( 9 )=TI LT*TT<  3 > 

DO  1 1=1,32 
J = ( I-l  1*2  + 1 
K= I T A(  I 1 

A(  I 1=AA(J1*TT(K1+AA(J+1  1*TT(K+2  1 
B(  I )=BB<J  }*TT(K  1 + BB(J  + 1 l*TT(K  + 2 ) 

K = 1TC(  I 1 

E(I  l=EE(Jl*TT(K)  + EE(J+l  l*TT{K  + 2) 

1 F(I  )=FF(J)*TT<K)+FF(J+1  1*TT(I(+2) 

DO  2 1=1,22 

J=(  I-l  1*2  + 1 
K=ITB(  I 1 

C(  I 1=CC(J1*TT(K1+CC(J+1  1*TT(K  + 21 

2 D(I  1=D0(J)*TT(K  1+DD(J  + 1 1»TT(K  + 21 

C CALCULATED  THE  FIELD  FR0P1  THE  POLYNOMIAL  EXPANSION 
6 EXPR  = EXP( -0 , 06*R2  1 

11  = 1 
JJ  = 1 
KK=  1 


K = 0 

10  BZ=BZ+(E(KK)+F(KK)»EXPR)»ZEYEXB 
KK=KK+1 

BX  = BX  + ( A(  1 IUB(  1 1 )*EXPR  )*ZEYEXB 
I 1 = 11  + 1 

IF(  I JK  .GT.  8 ) GO  TO  20 

BY  = BY  + (C(JJ  ) + 0(  JJ  )»EXPR  )»ZEYEXB*Y 

JJ= J J+l 

15  ZEYEXB=ZEYEXB*Z 

I JK=1 JK+1 
K = K*1 

!F(  IJK.LE.9. AN0.K.LE.5  ) GO  TO  10 
20  YEXB=YEXB*Y2 

30  XB=XB*X 

MO  BF(  1 )=BX*CON 

BF( 2 )=BY»CON 
BF( 3 )=BZ*CON 
RETURN 

C ERROR  EXIT 
50  WRITE(6,60)  XX 

60  FORnATfMH  X=  ,E10.3,9H  Y=  ,E10.3,MH  Z=  ,E10.3,T6H  IS  OUTSIDE  THE 
*VALID  REGION--POUER  SERIES  DIVERGES  BFIElD  IS  SET  TO  ZERO  ) 

GO  TO  MO 
END 


SUBROUHNE  SPHRCM 

C SPHRCH  IS  A nODIMED  VERSION  OF  J.C.  CAINiS  14  TERn  FAST  SPHRC 
C ROUTINE. 

C IT  HAS  BEEN  SHORTENED  TO  10  TERMS. 

C IT  HAS  A TRUNCATION  FOR  LARGE  R - THE  TRUNCATION  BETWEEN  TERMS 
C IS  SMOOTH. 

C IT  IS  SELF  contained  AND  SET  UP  TO  USE  GAUSS  NORMALIZED  COEFFICIENTS. 

C 

C THE  IGS75  GAUSS  NORMALIZED  COEFFICINTS  THRU  N=10  ARE  INCLUDED  WITH 
C THIS  DECK  ( BARR ACLOUGH, ET  AL  J.  GEOPV.  R E S . , 6 45 ( 1 9 75  ) . 

C THEY  ARE  valid  FROM  1945.  TO  1975.  THEY  MAY  BE  USED  FOR  LATER  PERIODS 
C WITH  CAUTION  (I.E.  THEY  HAVE  NOT  BEEN  VERIFIED) 

C ARGUMENT  TRASMITION  IS  VIA  COMMON. 

C INPUT  PARAMETERS 

C YEAR  IS  THE  YEAR,  IF  YEAR  CHANGES  THE  COEFFICIENTS  ARE  UPDATED. 

C ST  SINE  OF  THE  GEOGRAPHIC  CO-LATITUDE. 

C CT  COSINE  OF  THE  GEOGRAPHIC  CO-LATITUDE. 

C SPH  SINE  OF  THE  geographic  LONGITUDE. 

C CPH  COSINE  OF  THE  GEOGRAPHIC  LONGITUDE. 

C AQR  6371. 2/R,  WHERE  R IS  THE  GEOCENTRIC  DISTANCE  IN  KM  FROM  THE 

C CENTER  OF  THE  EARTH. 

C NMAX  MAXIMUM  NUMBER  OF  TERMS  TO  BE  USED  (MUST  BE  LESS  OR  EQUAL 

C TO  10).  THIS  ROUTINE  PRESETS  IT  TO  10 

C OUTPUT  parameters 

C BR  RADIAL  COMPONENT  OF  FIELD  IN  GAUSS. 

C BT  THETA  COMPONENT  (SOUTH  POINTING)  COMPONENT. 

C BP  PHI  COMPONENT  (EAST) 

01  MENS  ION  G(10,10),CONST(!0,10),FM(10),FN(10) 

DIMENSION  G0( 10,10),GT( 10,10),GTT(10,10) 

COMMON /BXYZ Cm/ YE AR , DAYYR,UT,XMLG,K00E,JSU,!N1TL,XMLT 

COMMON  /GCOM/  ST, CT,  SPH, CPH,  ADR,  BT,  BP,  BR,  NMAX  1-00005 

DIMENSION  CONA(  lO ) 

DATA  YRLAST,TF,TL,TZER0,NMAX/-12345.,1945.,1975.,1975.,10/ 
data  IFIRST/0/ 

DATA  CONA/0 :,  12. 0,6. 5, 4. 5, 3. 5, 2. 9, 2. 6, 2. 3, 2. 1,2.0/ 

C COEFFICIENTS  FOR  BASE  YEAR  (1975.) 

DATA  GO  /O., 30103. 6, 2860. 05, -3195. 5, -4142. 69, 1737. 23, -636. 69, -1917 
•.09,-553.01,-883.14,-5682.6,2016.5,-5213.3,6558.51,-4305.68,-3572. 
*54,-1321.33,1890.53,-341.86,-1274.03,3576.17,50.32,-1414.22,-2429. 
•72,-1736.64,-2015,83,-413.96,-66.61,145.81,-173.84,1009.8,-514.91, 

•179. 46, -656. 96,  844. 82, 300.26,  1935. 78, -274, 41, 521. 89, 946.  ,-1070. 27, 
•1040.11,-110.86,210.91,-157.15,349.42,4.91,79.03,368.96,-597.58,-2 
*49.08,-1140.49,759.11,185.03,-64.75,28.2,-8.84,-19.76,1.48,-20.21, 
•211.71,-1500  4,-773.11,219.91,18.38,-10.48,73.01,-41.17,16.48,3.48 
•,2716.97,715.33,92.15,-86.44,-151.27,52.8,8.35,3.82,-30.83,-4.52,- 
•328.45,779.54,-207.1,481.25,-34.53,-99.55,27.83,10.47,-3.07,-1.29, 

•2497.  11, -1705. 81, -406. 62,  1 74. 76,  141. 5, -168. 76, -91. 91,  .52, -.18, -.3/ 

C FIRST  TIME  derivatives  OF  THE  COEFFICIENTS 


DATA  GT  /O  , -26. 82, 37. 56, 9. 43, 3. 85, -1.58, -9.  38,  10.  19, -22.  12,0.  ,10. 
•1,-10. ,-.45, 32. 03, 12. 23, 10. 37, -IT. 58, 6. 74, -21. 44,0. ,4. 88, 16. 37, -4. 
*75, 9. 12, 15. 73, -9. 91, -34. 82, 14, 77, 2. 24,0. ,-22. Cl, -5. 4 8, 5 06,3.72,4. 
•43, 9. 84, -35, 27, -5. 32, -17. 4,0. ,-29. 77, -2. 66, -535, .5, 3. 39, 1.35, -.05 
•,-10, 37, 6.  15,0.  ,-9.35, -20.  37,  12. 52, -2. 97, -.79, -.91, -1.88, -3  46, 5. 4 
•9,0. ,6  . 24, 3. 29, -2, 29, 8. 51, -.41,-1  .34,  .27,-1  . 16,-4.26,0.  ,43.27,5  .79 
•, -1.02, -3. 83, 3. 89,,!,-. 78, .49, .68,0. ,12. 07, 19. 07, 12  01,8.56,-7.27, 
•3.43,1.4,-.31,.01,0.,0.,0,,0.,0.,0.,0.,0.,0.,0.,0./ 

C SECOND  TIME  DERIVATIVES  OF  THE  COEFFICIENTS 


r 


DATA  GTT/0.,-.35,.15,.35,0...51,-.5l,0.,0.,0.,.25,0.,0,,0.,.‘*7,.71 
•,-.38,0.,0.,0.,-.59,0.,-.07,.31,0.,.5H,-.9,0.,0.,0.,-.2,0.,0.,0.,. 
*15,.2M,~.55,0.,0.,0.,-.83.0.,0.,-.06,.06,0.,0.,0.,0..0.,.71,-.38,0 

*.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0 

*.,0.,0.,0.,0,/ 

C SET  UP  INITIAL  CONSTANTS  DURING  FIRST  CALL 
IF(  IFIRSr.NE.O  ) GO  TO  199 
IFIRST=l 
Fn(  1 ) = 0 
00  6 N=2, 10 
Fru  N )=N-l 


FN(  N )=N 
DO  6 n=l,N 

t>  C0NST(N,n)  = FL0AT((N-2  >«»2-(n-l  )»*2  )/FL0AT((2»N-3)*(  2*N-5  I ) 

C SET  UP  THE  COEFFICIENTS 

C IF  YEAR  MAS  CHANGED  BY  MORE  THAN  .lYEAR  UPDATE  THE  COEFFICINTS. 

199  IF(  ABS<  YRLAST-YEAR >. LT . 0 . 1 ) GO  TO  230 
IF(  YEAR. GE. TF. AND. YEAR. LE. TL  ) GO  TO  210 
WRITE!  6,200  ) VEAR 

200  FORMAT! 10H0*»WARNING, FT. 1 ,23H  IS  OUTSIDE  TIME  LIMITS) 

210  T=VEAR-TZERO 

DO  220  Nrl,NnAX 
DO  220  M=1,NMAX 

220  GIN,ni=GO!N,n)*T*!GT!N,M>+T*GTT!N,n)) 

YRL  AST  = YEAR 
C 

230  AR= AOR»AOR*AOR 
AR= AOR» AOR»AOR 
C2=G!2,2)*CPH*G!1,2)*SPH 
Bn  = -!AR*AR)*(G!2,I  >*CT»C2*ST) 

BT=AR*!C2*CT-G!2,1  )*ST) 

BP= AR*! G!l,2  )*CPH-G!2,2)*SPH) 

IF  (NMAX.LE.2  ) RETURN 
R=1 . / AOR 

IF! R . GT . CONA! 2 ) ) RETURN 
CON=0 . 

SP2=SPH 

CP2=CPH 

P21=CT 

P22=ST 

DP21=-ST 

DP22=CT 

N=3 

SP3=! SP2*SP2  )»CP2 

CP3  = ! CP2*-SP2  )•!  CP2-SP2  ) 

P31  = CT»P21-C0NST<  3, 1 ) 

P32=CT*P22 

P33=ST*P22 

DP31=-P32-P32 

DP32=CT*DP22-P33 

DP33=-DP31 

C2  = GI 3,2  )*CP2  + G( I , 3 )*SP2 
C3=G!3,3)»CP3*G!2,3)*SP3 
ARr AOR* AR 

XR:BR-FN( 3 )*AR»( G( 3, 1 I « P 3 1 * C2 • P 3 2 ♦ C 3 • P 3 3 ) 

XT  = BT*AR.!G!3,1)»DP31  + C2»DP32*C3'>0P33) 


1-00006 
1-00007 
1-00008 
1-00009 
1-00010 
1-000  I 1 


1-00012 

1-00013 

1-00019 

1-00015 

1-00016 

1-00017 

1-0001  9 
1-00020 
1-00021 
1-00022 
1-00023 
1-00029 
1-00025 
1-00026 
1-00027 
1-00028 
1-00029 
1-00030 
1-00031 


XP  = BP-AR*(Fn(2)*(G(3,2)*SP2-G(l,3)*CP2)»P21*FM(3)*(G(3,3>»SP3-G(?,l 

♦ 3 )*CP3  )*P22  ) 1 

BP=BP*ST  1 

XP=XP*ST  t 

IF(NMAX.LE.3  ) GO  TO  21 

IF( R. GT. C0NA( 31)  GO  TO  20 

SP*«  = SPH*CP3+CPH*SP3  1 

CPM=CPH*CP3-SPH»SP3  1 

P*<  1 = CT*P31 -C0NST(  R,  1 )*P21  I 

DP41  = CT»DP31-ST»P31-C0MST(  i*,  1 )*DP21  I 

P'<2=CT*P32-C0NST(  **,2  )*P22  I 

DPM2=CT*0P32-ST*P32-C0NST( ^ , 2 )*DP22  1 

P*f3  = CT*P33  1 

0PM3=CT*DP33-ST*P33  1 

PS'4  = ST*P33  1 

DP‘4H  = Fn(  H )*PR3  I 

C2  = G(  <4 , 2 )»CP2  + G(  1 , 4 )*SP2  1 

C3=G( 4 , 3 )*CP3  + G( 2, 4 )*SP3  1 

C4  = G(  4 , 4 )*CP4  + G(  3 , >*SP4  1 

AR=A0R*AR  I 

BR  = XR-F(\I(4)»AR*(G(4,1)*P41*C2*P42  + C3*P43*C4»P44)  1 

3T=XT  + AR*( G( 4,  1 )»DP4 1 + C2»0P4  2 + C3*DP4 3 + C4 *0P44  1 1 

BP  = XP-AR»(Fm2)»(G{4,2)»5P2-G(l,4)*CP2)*P42*Fn(3J*(G(4,3l*SP3-G(2,I 

♦ 4)»CP3)*P43  + Fn(4)*(G(4,4)*SP4-G(3,‘4)*CP4)*P44)  1 

IF( NMAX . LE.4  ) GO  TO  11 

IF( R. GT. CONA( 4 ))  GO  TO  10 


N=5 


SP5  = < SP3  + SP3  )*CP3 

CP5  = ( CP3  + SP3  )*( CP3-SP3  1 

P51  = CT»P4 1-C0NST( 5 , 1 )*P31 

DP51=CT»DP41-ST*P41-C0MST(5,1)»DP31 

P5  2=CT*P4  2-C0VST( 5,2  )*P32 

DP52=C-f  *DP4  2-ST»P42-C0NST(  5 ,2  )*0P32 

P53=CT*P43-C0NST(5,3)*P33 

0P5  3=CT*DP4  3-ST*P4  3-C0NST( 5, 3 )*0P33 

P54=CT*P44 

0P54=CT*0P44-ST*P44 

P55=ST*P44 

DP55  = FM(  5 )oP54 

C2=G(5,2  )*CP2  + G(  1,5  )*SP2 

C3  = G( 5 , 3 )*CP3  + G( 2, 5 )*SP3 

C4  = G(5,4  )»CP4  + G(3,5  )*SP4 

C5  = G(5,5  )»CP5*G<4,5  )*SP5 

Afl= A0R*AR 

XR  = BR-FN(5  )*AR*(G<5,1  )»P5  1+  C2*P5  2 + C3»P5  3+C4»P54  + C5*P55  ) 

XT=BT*AR»(G(5,1  )*DP51  + C2*DP52  + C3»DP53*C4»DP54*C5*DF55  1 

XP  = BP-AR.(Fn(2)*(G(5,2)»SP2-G(l,5>*CP2)»P52  + F(«U3).(G(5,3)»SP3-G(2, 

♦ 5 )»CP3  )»P5  3 + Fn(4  )*(G(5,4  )*SP4-G(3,5  )»CP4  )*P54  + FM(5  )»(G(5,5  )*5P5-G( 

♦ 4,5  )*CP5  )»P55  1 

IF( NMAX.LE.5  ) GO  TO  21 
1F( R. GT. C0NA( 5)1  GO  TO  20 


Nrfc 

SP6=SPH*CP5*CPH*5P5 

CP6=CPH»CP5-SPri*SP5 

P61  = CT»P51-C0NST(6,  1 )»P41 

DPfcl=CT•DP51-ST•P51-C0^ST(6, 1 >*0P41 


-00032 

-00033 

-00035 

-00035 


-0003T 
-00038 
-00039 
-00040 
-00041 
-00042 
-00043 
-00044 
-00045 
-00046 
-0004T 
-00048 
-00049 
-00050 
-0005  1 
-00052 
-00053 
-00054 


-0005T 

-00058 

-00059 

-00060 

-00061 

-00062 

-00063 

-00064 

-00065 

-00066 

-00067 

-00068 

-00069 

-00070 

-00071 

-00072 

-00073 

-00074 

-00075 

-00076 

-00077 

-00078 


-00081 

-00082 

-00083 

-00084 


P62=CT»P52-C0NST(6,2)»P42  1-00085 

0P6  2=CT*0P5  2-ST*P5  2-C0NST(  6,2)*0P‘«2  1-00086 

P6  3=CT*P5  3-C0NST(  6,3  )»P‘4  3 1-0008T 

DP6  3 = CT»0P5  3-ST*P5  3-C0NST( 6,3  )*DPM3  1-0008  8 

P6H  = CT*P5‘4-C0NST(  6,')  1-0008<> 

DP6‘*  = CT*DP5M-ST*P5‘*-C0NST(  6,‘l  l^OP***!  1-00090 

P65=CT*P55  1-00091 

DP65=CT*DP55-ST*P55  , 1-00092 

P66=ST*P55  1-00093 

DP66  = Pf1(  6 )*P65  1-00099 

C2=G( 6, 2 )*CP2*G( 1 , 6 )*SP2  1-00095 

C3  = G( 6, 3 )*CP3  + G( 2, 6 )*SP3  1-00096 

C9  = G( 6 ,9  )»CP9  + G( 3, 6 )*SP9  1-00097 

C5  = G(6,5  >*CP5  + G(9,6  )»SP5  1-00090 

C6  = G( 6,6  )*CP6  + G( 5,6  >»SP6  1-00099 

AR=AOR*ftR  1-00100 

BR=XR-FN(6)*AR*(G(6,1 )*P61+C2*P62+C3*P63+C9»P69+C5*P65+C6*P66  > 1-00101 

BT=XT+AR*(G(6,1)*DP61+C2*DP62+C3*DP63+C4*DP69+C5*DP65+C6*DP66)  1-00102 

BP=XP-AR*(FnU2)*(G(6,2)»SP2-G(l,6)»CP2)*P62*Fn(3>*(G(6,3)*SP3-G(2,l-00l03 
+ 6)*CP3>-»P63  + Fn(9)*(G(6,9)*SP9-G(3,6)*CP9)*P69  + Fn(5)*(Gt6,5)*SP5-G(I-00109 
♦ 9,6)»CP5)*P65t-Fn(6)*(G(6,6)*SP6-G(5,6)*CP6)*P66)  1-00105 

IF( NMAX . LE. 6 ) GO  TO  11 
IF( R . GT. CONA( 6 ) > GO  TO  10 
N=7 


SP7  = ( SP9  + SP9  )*CP9 

1-00108 

CP7  = ( CP9  + SP9  )*( CP9-SP9  ) 

1-00109 

P71  = CT*P61 -CONST! 7,1  )»P5l 

1-00110 

0P71=CT*DP61-ST*P61-C0N5T(7,1)*DP51 

1-0011  1 

P72=CT*P62-C0NST(7,2)*P52 

1-00112 

DP72=CT*DP62-ST*P62-C0NST(7,2)*DP52 

1-00113 

P73=CT»P63-C0NST(7,3)*P53 

1-00119 

DP 7 3= CT*DP6  3-ST*P6 3-CONST! 7, 3 >*DP5  3 

1-00115 

P79=CT*P69-C0NST(7,9)»P59 

1-001 16 

DP 79  = CT*DP69-ST»P69- CONST! 7,9  )*DP59 

1-00117 

P75=CT»P65-C0NST1 7, 5 )*P55 

1-00118 

0P75=CT*DP65-ST*P65-C0NSTI7,5)*DP55 

1-00119 

P76=CT*P66 

1-00120 

OP  76  = C 7* DP6 6 -ST*P6 6 

1-00121 

P77=ST*P66 

1-00122 

DP77  = FM( 7 )*P76 

1-00123 

C2=G!7,2)»CP2+G!1,7)»SP2 

1-00129 

C3=G!7,3)»CP3  + GI2,7)*SP3 

1-00125 

C9=G!7,9)*CP9+G(3,7)*SP9 

1-00126 

C5  = G! 7,5  )«CP5+G(9, 7 )»SP5 

1-00127 

C6=G(7,6)*CP6+G!5,7)*SP6 

1-00128 

C7  = G! 7, 7 )*CP7  + G! 6, 7 )*SP7 

1-00129 

AR= AOR»AR 

1-00130 

XR=BR-FN!7)*AR»!G!7,1)*P71+C2*P72*C3*P73+C9»P79+C5*P75+C6 

•P76+C7»P1-00131 

♦77)  1-00132 

XT=BT+AR.(G(7,l)»0P71+C2*DP72+C3*DP73+C9*DP79+C5*DP75+C6*DF76+C7*Dl-00133 
♦P77)  1-00139 

XP  = BP-AR»(Fn(2)*(G(7,2)*SP2-Gil,7)»CP2)*P72  + F,'0(3>*(G(7,3)»SP3-G(2;l-00135 
♦7)»CP3>»P73^Fn(9)»(G(7,9)»5P9-G(3,7)*CPM)*P79+Fn(5)*<Gf7,5)»SP5-G(l-00136 
♦ 9,7)»CP5)*P75  + Fm(6)*(G(7,6>»SP6-G(5,7)*CP6)*P76  + FPl(7)*(G(7,7)*5P7-l-00l37 
♦G( 6, 7 )*CP7 )*P77  ) 1-00138 

fF(  \nax.LE. T ) GO  ro  2i 
IF( R , GT. CONA( 7 ) ) GO  TO  20 


N=8 


SP8=SPH*CP7+CPH*SP7 

1-OOlHl 

CP8=CPH*CP7-5PH*SP7 

1-001M2 

P81  = CT*P71-C0NST( 8. 1 )*P61 

1-001H3 

DP81  = CT*DP71-ST*P71-C0NST(8,n*DP61 

1-001<<‘* 

P82  = CT*P72-C0rjST(8,2)*P62 

1-001*«5 

DP82=CT*DPT2-ST*P72-C0N5Tf8,2)*0P62 

1-001‘«6 

P83=CT»P73-C0NST(8,3)*P63 

l-OOl'^T 

DP83=CT*0P73-ST»P73-C0NST(8,3)*0P63 

l-00l'*8 

F8R  = CT»P7R-C0NST(8,R  )*P6H 

1 -001 

DP8R  = CT*DP7R-ST»P7R-CONST(8,R)*OP6‘* 

1-00150 

P85  = CT*P75-C0NST(8,5  )»P65 

1-00151 

0P85  = CT»DP75-ST*.'75-C0NST(8,5>*DP65 

1-00152 

P86  = CT»P76-C0NST(8,6  )*P66 

1-00153 

DP8  6 = CT*DP76-ST*P76-C0NST( 8,6  )*DP6  6 

1-0015** 

P87=CT»P77 

1-00155 

DP87=CT»0P77-ST*P77 

1-00156 

P88=ST»P77 

1-00157 

DP88  = Fn( 8 )*P87 

1-00158 

C2=G(8,2)*CP2*G(1,8)*SP2 

1-00159 

C3=G( 8, 3 )*CP3  + G( 2, 8 )*SP3 

1-00160 

Ca  = G(8,M  )*CPR  + G(3,8)*SPR 

1-00161 

C5  = G( 8,5  )*CP5+G(R,8  )*SP5 

1-00162 

C6=G(8,6)*CP6*G(5,8)*SP6 

1-00163 

C7=G(8, 7 )*CP7+G(6,8  )*SP7 

1-0016** 

C8  = G(8,8  )*CP8  + G(7,8  )*SP8 

1-00165 

AR=AOR*AR 

1-00166 

8R=XR-FN(8)*AR*{G(8,1)*P81+C2*P82*C3*P83*CR 

• P8***C5'*P85*C6*P86*CT*Pl-0016r 

♦ 87t-C8*P88  ) 1-00168 

BT=XT+AR*( G(8,1)*OP81+C2*DP82+C3*DP83*Cm*DP8M+C5*OP85*C6*DP86*CT*D1-00169 

+ P87  + C8*DP88  ) 1-00170 

BP=XP-AR*(FnU2)*(G(8,2)*SP2-G(l,8)»  CP  2)*?82  + Fn(3)*(G(8,3)»SP3-G(2, 1-00171 

♦ 8)*CP3)*P83+Fn(‘4)*(G(8,R)*SPR-G(3,8)*CPA)*P8M*Fn(5)*(G(8,5)»SP5-G(l-00172 
+R,8)*CP5)*P85+Fn(6)*(G(8,6)*SP6-G(5,8)*CP6)*P86+Fn(7)*(G<8,7)*SP7-l-00173 

♦ G(6,8)*CP7)*P87  + Fn(8)*(G(8,8)*SP8-G(7,8)*CP8)*P88)  1-0017** 

IF( NHAX . LE . 8 ) GO  TO  11 


IF(R.GT.C0NA(8  ) ) GO  TO  10 
N=9 


SP9  = ( SP5>SP5  )«CP5 

1-00177 

CP9  = ( CP5  + SP5  >*( CP5-SP5  ) 

1-00178 

P91  = CT*P81-C0NST( 9, 1 )*P71 

1-00179 

DP91=CT»DP81-ST*P81-C0NST(9,1)*0P71 

1-00180 

P92=CT*P82-C0NST(9,2)»P72 

1-00181 

DP92=CT*0P82-ST*P82-C0NST(9,2)**DP72 

1-00182 

P93=CT»P83-C0NST(9,3)»P73 

1-00183 

DP93  = CT*DP83-ST*P83-C0NST(9,3)'»DP73 

1-0018** 

P9**  = CT*P8‘*-C0NST(  9,**  )*P7** 

1-00185 

DP9**  = CT*DP8‘*-ST*P8‘*-C0NST(9,**)*DP7** 

1-00186 

P95  = CT»P85-C0NST( 9,5  )»P75 

1-00187 

DP95  = CT»DP85-ST*P85-C0NST( 9,5  )»DP75 

1-00188 

P96  = CT*P86-C0NST( 9,6  )*P76 

1-00189 

DP96=CT»0P86-ST*P86-C0NST(9,6)*DP76 

1-00190 

P97=CT»P87-C0NST<9,7)*P77 

1-00191 

DP97=CT*DP8 7-Sr*P87-C0NST{ 9, 7 )*DP77 

1-00192 

P98  = CT**P88 

1-00193 

DP98=CT»DP88-ST*P88 

1-0019** 

P99=ST*P88 

1-00195 

t 


0P99  = Fn( 9 )*P98  1-00196 

C2  = G(9,2  )»CP2*G( 1,9  )*SP2  1-00197 

C3=G( 9, 3 )*CP3*G( 2, 9 )*5P3  1-00198 

C9  = G( 9, 9 )*CP9*G(  3, 9 )*SP9  1-00199 

C5  = G( 9,5  )*CP5  + G(9,9 »*SP5  1-00200 

C6  = G( 9, 6 )*CP6  + G( 5 ,9  )»SP6  1-00201 

Cr  = G( 9, 7 )*CP7»G( 6,9  )»SPT  1-00202 

Ce=G( 9,8  }*CP8*Gl 7,9>*SP8  1-00203 

C9  = G( 9,9  )*CP9»G( 8,9  )»SP9  1-00209 

ftR=ftOR*AR  1-00205 

XR=BR-FN(9)*AR*(G(9,1)»P91*C2*P92*C3»P93+C9»P99*CS»P95*C6»P96*C7*P1-00206 

♦ 97t-C8*P98*C9»P99  ) 1-00207 

XT=BT*AR*(G(9,l)»DP91»C2»0P92*C3*0P93*C9»0P99*C5»0P95*C6»0P96*C7»0l-00208 

♦ P97  + C8*DP98*C9*DP99  ) 1-00209 

XP=BP-AR*(  Fm  2>*(G(9,2)*SP2-G(1,9)»  CP  2)»P92*FH(3)*(G(9,3)»SP3-G(2, 1-00210 

♦ 9)*CP3)»P93  + Fn(9)»(G(9,9)»SP9-G(3,9>*CP9>*P99*Fn{5)»<G<9,5)»SP5-G(l-00211 

♦ 9,9)*CP5)*P95»Fn{6)*(G(9,6)*SP6-Gf5,9)»CP6)»P96*Fm7>*«G(9,7)»SP7-l-00212 

♦ G(6,9)*CP7)»P97  + FM(8)*(G<9,8)»SP8-G(7,9)»CP8)»P98*F(V9)»1G<9,9)*SP1-00213 

♦ 9-G(8,9  )*CP9  )*P99  ) 1-00219 

1F( NnAX.LE.9  ) GO  TO  21 

IF( R. GT. C0NA( 9 ) ) GO  TO  20 
N=10 

SP10=SPH*CP9*CPH*SP9  1-00217 

CP10=CPH»CP9-SPH*SP9  1-00218 

PI  01  = CT*P91 -CONST!  10,1  )*P81  1-00219 

OPl 01=CT»DP91-ST»P91 -CONST! I0,1)»DP81  1-00220 

P102  = CT*P92-C0NST( 10,2  )*P82  1-002  21 

Dp102=CT*DP92-ST*P92-C0NST(10,2)*0P82  1-00222 

P103=CT*P93-C0NST(10,3>*P83  1-00223 

DPI 03=CT*OP93-ST»P93-CONST( 10,3)*DP83  1-00229 

P109=CT*P99-C0NST(10,9)*P89  1-00225 

DP109  = CT»DP99-ST*P99-C0NST(  10,9  )»0P89  1-00226 

P105  = CT*P95-C0NST(  10,5  )*P85  1-00227 

DP105  = CT*0P95-ST*P95-C0NST(  10,5  >*0P85  1-00228 

P106=CT»P96-C0NST(10,6)*P86  1-00229 

0P106=CT*DP96-ST*P96-C0NSTtl0,6)*0P86  1-00230 

P107=CT*P97-C0NST( 10,7)*P87  1-00231 

DP107=CT*DP97-ST»P97-CONST(10,7 )*0P87  1-00232 

PI  08  = CT*P9 8-CONST!  10,8)»P88  1-00233 

D?108=CT*DP98-ST*P98-C0NST!10,8)»DP88  1-00239 

P109=CT»P99  1-00235 

DP109=CT»DP99-ST*P99  1-00238 

P1010=ST*P99  1-00237 

DPIOIO=Fffl( 10  )*P109  1-00238 

C2  = G! 10,2  )*CP2*G! 1,10  )*SP2  1-00239 

C3=G!10,3)*CP3*GI2,10  )*SP3  1-00290 

C9=G!10,9)*CP9+G!3,10)»SP9  1-00291 

C5=G!10,5)»CP5*G!9,10)*SP5  1-00292 

C6  = G!  10,6  )»CP6  + G! 5,10  )*SP6  1-0029  3 

C7=G! 10, 7 )»CP7  + G! 6, 10  )*SP7  1-00299 

C8  = G! 10, 8 )»CP8*G! 7,10  )*SP8  1-0029  5 

C9=G!10,9)»CP9+G!8,10)*SP9  1-00296 

ClO  = G! 10, 10  )*CP10  + G!9, 10  )*SP10  1-00297 

AR=AOR*AR  1-00298 

BR=Xfl-FN!10)*AR»!G!10,l  )*P101*C2«P102  + C3*P103+C9»P 109* C5»P105*C6*P 1-00299 

♦106+C7*P10T*C8*P108+C9»P1C9*C10»P1010)  1-00250 

8T=XT* AR»! G! 1 0,1)* DPI 01  + C2* DP 102*C3»DP103+C9» DPI  09* C5» DPI  05*C6*DP1  1-0025  1 


73 


c 


♦ 06*C7»DP)!)7*C8*DP]09  + C'}*DP109*ClO»OPt010)  1-002S2 

BP=XP-AR*(Fn(2)»(G(10,2)*SP2-G(l,10)»CP2)»Pl02*Fn(3)»(G{lO,3)*SP3-l-00253 

♦ G(  2,  10  )*CP3  )»P103-»-Fn(  9 )•(  G(  10, 9 )»SP'i-0(  3,  10  )»CP9  )*Pi09*Fn(  5 )*(  G(  101-0025** 
+ ,5  )*SP5-G(**,10  )*CP5  )»P105*Fm6  )♦(  G(  10, fc  )*SP6-G(5,1J  )»CPfe  )*P106*Fn<  1-00255 

♦ 7 )*(G(  10,7  )»SP7-G(6,  10  )»CP7  )»P107*-F/n(8  )»(G(  10,8  >*SP8-G(7,10  )»CP8  )»l-00256 

♦ PlOe  + FCKR  )*(G<  10,9  )*SP9-G(8,  10  )<>CP9  )»P109*Fn(  10  )*(  G(  10,10  )*SP10-G(  1-00257 

+ 9, 10  )*CP10  )*P10  10  ) 1-00258 

I BP  = BP/ST  l-00**7‘* 

IF(  NPIAX  . LE  . 1 0 ) RETURN 

WRITE  ( 8,2  ) NHAX  1-00**76 

STOP  l-00<*77 


2 F0RnAT(57H0  ERROR,  THIS  SPHRC  ONLY  FOR  NnAX  = <'l0,  CALL  WAS  FOR  NMAX 
•=,  15  ) 

C MAKE  A SnOOTH  FIT  BETWEEN  TRUNCATED  TERPIS. 

10  C0N=(R-C0NA(N))/(C0NA(N-1  )-C0NA(N)) 

1 1 BR=BR^(  XR-BR  )»CON 
BT=BT*(  XT-BT  )*CON 

BP=( BP*(  XP-BP  )»CON  )/ST 
RETURN 

20  C0N=(R-C0NA(N))/iC0NA(N-1  )-C0NA(N)) 

21  BR=XR  + ( BR-XR  )*C0N 


BT=XT+( 8”-XT  )»CON 


BP=(XP+(BP-XP)*C0N)/ST 

RETURN 

END 


1-00*480 


subroutine  ANGLE(TILT,S1NPHE\,C0SPHE) 


c 

C THIS  ROUTINE  CALCULATES  THE  ANGLE  BETWEEN  THE  DIPOLE  AXIS  AND  THE 
C SUN-EARTH  LINE  AS  WELL  AS  THE  ROTATION  SINES  ANO  COSINES  TO  CONVERT 
C PRON  GEOMAGNETIC  TO  SOLAR  MAGNETIC  COORDINATES. 

C INPUT  PARAMETERS  ARE  VIA  COMMO N / B X V Z C M / 

C OAVYR  IS  THE  DAY  OF  YEAR  (1-366.).  IT  MUST  BE  A WHOLE  NUMBER. 

C DAY  1 IS  JANUARY  1 

C UT  IS  THE  UNIVERSAL  TIME  IN  HOURS  (0.-2M.) 

C OUTPUT  PARAMERS  ARE  VIA  THE  PARAMETER  LIST  IN  THE  ARGUMENTS 
C TILT  IS  90-P5I.  WHERE  PSI  IS  THE  ANGLE  BETWEEN  THE  DIPOLE  AXIS 
C AND  THE  SOLAR  DIRECTION. 

C S I NPHE  S I N( BET  A ) 

C COSPHE  COS(BETA)  BE'^A  IS  THE  ROTATION  ANGLE  FROM  GEOMAGNE'^IC 
C TO  solar  magnetic  COORDINATES 

COMMON/BXYZCM/YEAfi^  OAVYR,  UT,XMLG,K0DE,JSU,INITL,XMLT 
DATA  IFIRST/OV 

C THE  FIRST  TIME  THROUGH  THE  SUBROUTINE  SET  UP  FIXED  CONSTANTS. 

IF(  IFIRST . NE. 0 ) GO  TO  10 
IFIRST=1 

P 12= ATAN2( 0 . , -1  . )/2  . 

C0N  = 90 . /P  12 

SALF  = SIN(  1 1 . r/CON  ) 

calf=cos( 1 1 . r/CON  ) 

SG AM=S I N( 23 . 5 /CON  I 
CGAM=COS( 23 . 5 /CON  ) 

S ASG=S ALF*SG AM 
SACG=SALF*CGAM 
CASG=C ALF*SG am 
CACG=CALF*CGAM 
W=Pl2/6 . *(  1 . +1  . /365  .256  > 

C MAIN  ENTRY  SET  UP  SINES  AND  COSINES 
10  WT  = ( UT-1 6 . 6 + ( DAYYR-1 72  . )*2H.  )*W 

CWT  = -WT/365 .256 
SSMLT  = SI N( WT  ) 

CSMLT  = COS( WT  I 
SBWT  = SIN( CWT  ) 

CBWT  = COS<  CWT  ) 

SMLSWT=SSMLT*SBWT 

5MLCWT=SSMLT*CBWT 

cmlswt=csmlt*sbwt 

CMLCWT=CSMLT*CBWT 

C COMPONENTS  OF  THE  DIPOLE  AXIS  ARE  IN  ECCLIPTIC  COORDINATES 
ZX  = CASG*CBWT*S ACG»CMLCUT-S  ALF»SMLSWT 
ZVrCA5G*SBUT+SACG*CMLSWT+5ALF»SMLCWT 
2Z=CACG-SASG*CSMLT 
PSI=ACOS(ZX) 

0SP=1  . /(  SINCPSI  ) ) 

TILT=CON«( PI2-PSI  ) 

C COMPONENTS  OF  GEOMAGNETIC  X AXIS  IN  ECCLIPTIC  COORDINATES 
XX=CACG*CMLCUT-S ASG*C BUT- CALF. SMLSUT 
X Y=CACG»CMLSWT-S  flSG.SBWT.C ALF.SMLCUT 
XZ=-CASG*CSMLT-SACG 
SINPHE=(XY.ZZ-XZ.ZV).OSP 

COSPHE=( XX.(ZZ.ZZ.ZYZY)-ZX.(XV.2V*XZ.ZZ)).05P 

RETURN 

END 


nor) 


r 


subroutine  INTERPI BB, CC, D, E, J ) 

c 

C INTERPOLATION  ROUTINE 

C GIVEN  A SET  OF  THREE  POINTS  BB(3>  AND  CCI 3 ),  WITH  CC  THE  INDEPENDENT 
C VARIABLES  X AND  BB  THF  DEPENDENT  VARIABLE  V,  INTEflP  FINDS  THE  SOLUTION 
C TO  THE  THREE  LINEAR  EQUATION  EXPRESSING  V AS  A SECOND  ORDER 
C POLYNOnlAL  OF  X. 

C IF  J LESS  THAN  ZERO,  SET  0 TO  THE  rilNinun  VALUE  OF  THE  DEPENDENT 
C VARIABLE 

C IF  J GREATER  THAN  ZERO  SET  D TO  THE  VALUE  OF  THE  INDEPENDENT 
C variable  such  that  the  DEPENDENT  VARIABLE  V EQUALS  E. 

SINCE  THERE  ARE  TWO  ROOTS  TO  THE  EQUATION  CHOOSE  THE  ONE 
CLOSEST  TO  CC(K),  WHERE  K=3  IF  J=l,  K=2  IF  J=H,  K=1  IF  Jr2. 

DinENSION  BB(3),CC(3) 

Yl  = ALOG( BB( 1 ) ) 

Y2  = AL0G( BB( 2 ) > 

V3=AL0G( BB( 3 ) ) 

X2  = CC( 2 )-CC( 1 ) 

X3=CC( 3 )-CC( 1 ) 

DD=( X3-X2  )*X2*X3 
A=(X3*(Y1-Y2  ) + X2*(Y3-V1  ))/DD 
B = <X3**2*(Y2-Vl  >-X2**2*(Y3-Y1  ))/DD 
IF( J.LT.O  )GO  TO  100 
C=Y1-AL0G<  0 > 

0IS=B**2-*«  . *A*C 
IF( OIS . LE. 0 . ISTOP 
OIS  = SQRT( OIS  ) 

S.a=(-B  + 0IS)/(2.*A)+CC(1  ) 

SB=(-B-0IS)/(2.*A)-t-CC(I  ) 

E = SA 
K = 3 

IF( J.EO.M  ) K = 2 
IF( J . EQ . 2 ) K=1 

IF(ABS(SB-CC(K)).LT.ABS(SA-CC(K)))E=SB 
RETURN 
100  X=-B/(2.*A) 

E=EXP( A*X**2*B*X  + Y1  ) 

RETURN 

END 
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SUBROUTINE  DAXFI X( BMftX, BnIN  ) 

C 

C MAXFIX  nODIFIES  BMAX  SUCH  THAT  THE  COMPUTED  L BETTER  F[TS  THE 
C DRIFT  SHELLS  OF  ISOTROPIC  PARTICLE  DISTRIBUTIONS. 

C 

COnnON/BX YZCn/ YEAR, DAY VR , UT, XnLG.KODE.JSU.rNITL.xnLT 
C IF  JSU  LESS  than  OR  EQUAL  TO  ZERO,  AVERAGE  L SHELLS  ARE  NOT  WANTED. 
IF( J5U. LE . 0 IRETURN 

C IF  BMAX  IS  greater  THAN  3*BnlN  NO  MODIFICATION  IS  REOUIRED,  DRIFT 
C SHELLS  ARE  ALREADY  CLOSE  TO  AVERAGE  AND  NO  FURTHER  CHANGE  WILL  OCCUR. 
IF( BMAX . GT . 3 . »BMIN  ) RETURN 

C CALCULATE  PRELIMINARY  VALUE  OF  L — IF  L LESS  THAfg  3,  AVERAGE 
C SHELL  DOES  NOT  CHANGE. 

EL  = (0.311653/BnlN)**(l./3.  ) 

IF( EL.LT.3.  > RETURN 
CON=l  . 3 

IF(EL.LT.5.>  CONrl . +0. 15*( EL-3.  ) 

BrnAX=AMINl(3.»8MIN,C0N*BnAX) 

RETURN 

END 


\i 


w 


iH 


15 


6 

s 

1 0 
1 1 
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SUBROUTINE  CARHEL  (8,XI,VL>  CARniOOO 

COMPUTE  L CarmlOOI 

IF( xl-1  ,0E-3fc  MR,  IM,  15  CarmlOA? 

VL=(0.311653/B)»»(l./3.)  CARMLOB? 

RETURN  CARMLOCz 

XX  = 3 . 0*AU0G( X I ) 

XX=XX*AL0G(B/0.311653) 

IF(  XX*22 . )1 , 1 , 8 CARHLO03 

IE(  XX  + 3 . 12 , 2, 9 CARMLOOR 

IF(XX-3.  >3,3,10  CAPnL005 

I E(  X X-1  1 . 7 )R , R , n CARML006 

I  E(  X X-2  3 . >5,5,6  CARML007 

GGr . 333338*X X* . 30062 1 02  CARMloOB 

GO  TO  T CARHL009 

2 GG=((((((((-8.1537735E-1R*XX»8.3232531E-13>»XX*1.0066362E-9)»XX*  CarhloIO 

18.  10R8663C-8  >*XX  + 3.291635RE-6  >*XX  + 8.2711096E-5  >*XX*1  .371R667E-3  >♦  CAR ML 01 1 
2XX*.0150172R5>*XX+.R3R326R2>»XX+. 62337691  CARML012 

GO  TO  7 CARML0I3 

3 GG=((((((((2.60R7o23E-10*XX*2.3028767E-9)*XX-2.1997983E-8>»XX-  CARMLOIR 

15.39776R2E-7  >»XX-3.3R08822E-6  >»XX^3.8379917E-5  >*XX  + 1 . 178R23RE-3  >♦  CARML015 
2XX+1.RR92RR1E-2>»XX+.R3352788>»XX+.62286RR  CARML0I6 

GO  TO  7 CARML017 

9 GG  = ((((((((6.3271665E-10*XX-3. 958306 E-8)*XX  + 9.97661R8E-071»XX-  CARMLO I 8 

11  2531932E-5  >*XX  + 7.9R51313E-5  >»XX-3.2077032E-R  >»XX*2.  1680398E-3  >»  CARML0I9 
2XX»1.2S17956E-2>»XXi-.R3510529)»XX+. 6222355  CARML020 

GO  TO  7 CARML021 

5 GG=(((((2.8212095E-8*XX-3.80R9276E-6>»XX42.17022RE-R>»XX-6.7310339CARML022 

1E-3>*XX*.1203822R>*XX-.18R61796>*XX  + 2. 0007187  CAR ML  02 3 

GO  TO  7 CARML02R 

6 GG  = X X-3 . 0R6068  1 CARML025 

7 VL=(  ( ( 1 . 0 + EXP  ( GG > >»0 . 31  1 653  )/B  )»*(  1 . /3  . > CARML026 

END  COMPUTE  L CARML027 

RETURN  CARML028 

END  CARML029 
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5.4  Sample  Program  Output 
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Section  6.0 


APPENDIX  B - COORDINATE  TR;'\NSF0Rf1ATI0NS  AND  RELATED 
ANALYTIC  DERIVATIONS 

The  magnetospheric  magnetic  field  model  described  in  Section  2 is  represented 
in  solar-magnetic  coordinates.  To  use  the  model  it  is  frequently  necessary  to 
express  the  components  of  the  vector  field  in  geographic  or  geomagnetic  coor- 
dinates. Transformations  betv/een  these  coordinate  systems  are  described  below. 
Normally  such  transformations  are  time  consuming  and  non  analytic  because  of  the 
earth's  elliptical  orbit  around  the  sun.  Here  however  the  approximation  of  a 
circular  orbit  is  made.  This  permits  the  transformations  to  be  analytic  and  to 
require  a minimum  of  computer  time.  The  transformations  are  used  in  subroutines 
ANGLE  and  BI-INEXT  described  in  Appendix  A. 

6. 1 Coordinate  Transformations 

Define  a as  the  angle  betv/een  the  rotation  and  geomagnetic  dipole  axes,  y as  the 
angle  betv/een  the  rotation  axis  and  the  perpendicular  to  the  ecliptic  plane.  Tv/o 
angular  velocities  are  also  defined,  n as  the  velocity  of  the  earth  around  the 
sun  and  to  as  the  earth's  angular  velocity  about  its  rotation  axis.  Also  the  time, 
t,  is  measured  from  magnetic  noon  at  summer  solstice  when  the  rotation  and  dipole 
axes  and  the  sun-earth  line  are  coplanar. 

It  is  desired  to  add  the  external  field  (given  in  solar  magnetic  coordinates) 
vectorial ly  to  the  internal  field  represented  in  geographic  or  geomagnetic 
coordinates.  To  do  so  the  following  three  coordinate  systems  are  defined. 
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Ill  1 .1 

(1)  the  X Y Z ' ecliptic  coordinates  with  X always  toward  the  sun,  Z 
perpendicular  to  the  ecliptic  plane. 

(2)  the  XG,  YG,  ZG  coordinates  - with  ZG  coincident  with  the  geographic 
rotation  axis.  Tine  is  measured  such  that  at  t = 0 to  the  YG  axis  is 
in  the  ecliptic  XZ  plane  as  is  ZG.  The  angle  y between  ZG  and  Z is 
23.5“. 

(3)  the  Xti,  YH,  ZM  coordinates  - with  ZH  coincident  with  the  geomagnetic 
dipole  axis.  The  angle  a between  ZM  and  ZG  is  11.7°.  At  t = 0,  XM 
and  YH  are  in  the  XG,  ZG  plane. 

I 

Note  that  Y ,YG,  and  YM  are  not  usually  coplanar. 

The  external  field  is  given  in  4 variables,  the  position  X,Y,Z,  and  the  "tilt" 

I 

angle,  u,  between  the  X (ecliptic  axis  - toward  the  sun)  and  the  solar 
magnetic  Y axis  - defined  below. 

Solar  magnetic  coordinates  are  defined  in  terms  of  the  magnetic  dipole  axis  and 

t 

the  sun  earth  line  ( the  X axis  in  solar  ecliptic  coordinates).  The  Y^^^^  axis 

I 

is  in  the  plane  formed  by  the  X and  ZH  (=  Z^^^)  axes.  The  Y^j^  axis  then  forms 
the  orthogonal  set. 

Inputs  to  the  transformations  must  include  the  position  in  geographic  coordinates 
where  B is  to  be  found,  the  universal  time  of  day,  t,  and  the  day  of  the  year,  T. 

In  order  to  agree  with  these  coordinate  definitions  t and  T must  be  "corrected"  - 
to  magnetic  noon  summer  solstice.  T and  t as  discussed  here  are  "corrected"  - i.e., 
measured  from  magnetic  moon  at  summer  solstice. 


Normally  these  transformations  take  much  computer  time  and  involve  the  equation 
of  time  which  is  not  in  closed  form.  However,  great  computational  effort 
is  avoided  if  it  is  assumed  simply  that  the  earth  is  in  circular  orbit  around 
the  sun.  Then  all  transformations  are  analytic  and  require  little  computer  time. 

Given  X,Y,Z  in  geographic  coordinates,  and  t and  T (T  is  an  integer  variable), 
the  position  must  be  transformed  to  solar  magnetic  coordinates  and  the  tilt 
angle  found. 

X = r sin  6 cos  <p 
Y = r sin  6 sin  i}) 

Z = r cos  0 

where  0,<>  are  geographic  latitude  and  east  longitude  (measured  from  Greenwich). 


The  transformation  then  proceeds  as 
Rotate  to  the  dipole  longitude. 

XA  = X cos  6-  Y sin  6 
YA  = Y cos  6+  X sin  6 
ZA  E Z 

Rotate  to  the  dipole  axis  about  YA  axi 
XB  = XA  cos  a-  ZA  sin  a 
YB  = YA 


ows. 


XA 


ZB  = ZA  cos  a+  XA  sin  a 


r 


XB,  YC,  ZB  then  give  the  position  in  georr:etric  coordinates  and  the 
geomagnetic  longitude  of  the  point  is 


= tan-’(f). 


A final  rotation  about  the  magnetic  axis  is  needed  to  transform  to  solar 
magnetic  coordinates.  To  do  this  the  tilt  angle,  p,  must  be  computed.  This  is 

A\ 

best  done  by  expressing  H,  the  unit  vector  along  the  magnetic  dipole  axis,  in 
ecliptic  coordinates 


using  the  coordinates  defined  above. 


to  sun 


In  terms  of  unit  vectors  (their  direction  indicated  by  subsc-ipt),  here  the 
geomagnetic  coordinates  are  to  be  expressed  in  terms  of  ecliptic  coordinates.  Thus 
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5yM  “ ^Z6  " ^Z6  ^ 

^^ZM  “ ® ^XG  w t sin  a + J^q  sin  u t sin  a 

/k  A ^ ^ 

*^XM  * ^XG  ^ ° * ^YG  sin  u t cos  a - k^Q  sin  d 

/M  A|  I 

Then,  in  ecliptic  coordinates  (i,  j,  k ) 

^XG  “ - k'  sin  Y + cos  J1  t cos  y + 3 sin  J2  t cos  y 

jyg  = 3'  cos  fi  t - i'  sin  51  t 

k,>  = k*  cos  Y + 3'  sin  Y cos  n t + j sin  y sin  51  t 

ZG 

Then,  the  angle,  c.  between  the  earth  sun  line  and  the  dipole  axis  (the  complement 
of  the  tilt  angle,  y)  is  given  by 

cos  c » • 3*  which  by  inspection  is: 

cos  5 = cos  a sin  Y cos  51  t 

+ cos  (i)  t sin  o cos  51  t cos  y 

- sin  (1)  t sin  a sin  51  t 

Thus  c = cos"^  (cos  a sin  Y cos  51  t 

+ sin  a (cos  w t cos  51  t cos  y - sin  u t sin  fl  t)) 

and  finally  v * j ■ C* 

It  remains  to  find  the  rotation  angle  5 between  the  magnetic  prime  meridian  at 

A A A I 

the  equator,  i^^g,  and  the  solar  magnetic  X axis  (iQ|j|)-  P'"0"»  y ^ 

construct  the  (solar  magnetic  axis).  (M^j^j  = M^j^) 


or 


^ All 

- 

sin  c 


1 

sin  5 


1 

sin  c 


j 

0 


The  angle,  <b,  between  and  the  axis  (in  geomagnetic  coordinates  is  given 

by 

JsH  • Jyh  = 5i1t  «Z  - «y  ^ ♦ 


where 

My  = cos  a sin  y sin  « t + cos  w t sin  a sin  n t cosy+  sin  ut  sin  a cos  nt 

= cos  a cos  y - cos  w t sin  a sin  y 

6My  = cos  u t cos  £1  t - sin  u t sin  £1  t cos  y 
GM^  = sin  0)  t sin  y 

To  determine  (f  (the  angle  between  the  solar-magnetic  X axis  and  the  geomagnetic 

X axis)  over  the  full  range  of  angles  (0  to  2ir  radius)  it  is  necessary  on  the 

* m 

computer  to  use  the  tan'^  function  in  place  of  the  cos"  function  described  above. 
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A ^ 

Thus  (described  above  as  ly,.)  must  be  projected  onto  the  solar  magnetic 
GM 

A A 

and  axis.  Then, 

A A 

1/.U  “ leu  (cos  u t cos  a cos  n t cos  y - sin  u t cos  a sin  n t 
un  oil 

- sin  a sin  y cos  n t)  + 

jjj^  (cos  u)  t cos  a sin  n t cos  y + sin  u t cos  a cos  n t 

- sin  o sin  Y sin  t)  + 

A 

*^SM  u t cos  a sin  y - sin  a cosy) 

A A 

Then  Igj^  • l^^jj  * cos  w t cos  a cos  n t cos  y 

- sin  0)  t cos  o sin  t - sin  a sin  y cos  fl  t 

A A 

^GM  * ^SM  “ COS  u)  t cos  o sin  n t cos  y 

+ sin  u t cos  o cos  J1  t - sin  ot  sin  y sin  U t 


The  angle  41  Is  then  given  by 


♦ 


1 • 1 
GM  '’SM 

^6M  * ^SM 


The  transformation  betv/een  geomagnetic  and  solar  magnetic  coordinates  is  then 
simply  a rotation  through  the  angle  ^ (as  determined  above). 


Z3  £T  2^^ 


Given  the  magnetic  field  in  solar  magnetic  coordinates. 


6.2  Time  Derivation  of  the  Tilt  Anole 


Since  the  induced  electric  field,  E,,  resulting  from  the  daily  variation  in  the 
tilt  angle,  y,  is  given  in  terms  of  the  magnetic  vector  potential.  A,  (Ej  = - 
it  is  necessary  to  compute  the  time  variation  in  y.  This  is  because  "A  is  given  as 
a function  of  y and  not  t.  Thus 


y ■ ^ ^ » <t>  • cos“^  (ZX)  where 

ZX  cos  a sin  y cos  n t + cos  u t sin  a cos  a t cos  y 

- sin  (1)  t sin  a sin  fi  t. 


It  follows  that 

It  “ ■ It  * s^'ice  generally 

dY 

d (cos-'  (Y)) . 
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